{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TensorFlow version: 2.0.0\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "import datetime, os\n",
    "#hide tf logs\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'  # or any {'0', '1', '2'},\n",
    "#0 (default) shows all, 1 to filter out INFO logs, 2 to additionally filter out WARNING logs, and 3 to additionally filter out ERROR logs\n",
    "import scipy.optimize\n",
    "import scipy.io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker\n",
    "import time\n",
    "from pyDOE import lhs         #Latin Hypercube Sampling\n",
    "import pandas as pd\n",
    "import seaborn as sns \n",
    "import codecs, json\n",
    "from scipy.optimize import minpack2\n",
    "\n",
    "# generates same random numbers each time\n",
    "np.random.seed(1234)\n",
    "tf.random.set_seed(1234)\n",
    "\n",
    "print(\"TensorFlow version: {}\".format(tf.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Data Prep*\n",
    "\n",
    "Training and Testing data is prepared from the solution file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_1 = np.linspace(-1,1,256)  # 256 points between -1 and 1 [256x1]\n",
    "x_2 = np.linspace(1,-1,256)  # 256 points between 1 and -1 [256x1]\n",
    "\n",
    "X, Y = np.meshgrid(x_1,x_2) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Data\n",
    "\n",
    "We prepare the test data to compare against the solution produced by the PINN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_u_test = np.hstack((X.flatten(order='F')[:,None], Y.flatten(order='F')[:,None]))\n",
    "\n",
    "# Domain bounds\n",
    "lb = np.array([-1, -1]) #lower bound\n",
    "ub = np.array([1, 1])  #upper bound\n",
    "\n",
    "a_1 = 1 \n",
    "a_2 = 1\n",
    "\n",
    "usol = np.sin(a_1 * np.pi * X) * np.sin(a_2 * np.pi * Y) #solution chosen for convinience  \n",
    "\n",
    "u = usol.flatten('F')[:,None] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trainingdata(N_u,N_f):\n",
    "    \n",
    "    leftedge_x = np.hstack((X[:,0][:,None], Y[:,0][:,None]))\n",
    "    leftedge_u = usol[:,0][:,None]\n",
    "    \n",
    "    rightedge_x = np.hstack((X[:,-1][:,None], Y[:,-1][:,None]))\n",
    "    rightedge_u = usol[:,-1][:,None]\n",
    "    \n",
    "    topedge_x = np.hstack((X[0,:][:,None], Y[0,:][:,None]))\n",
    "    topedge_u = usol[0,:][:,None]\n",
    "    \n",
    "    bottomedge_x = np.hstack((X[-1,:][:,None], Y[-1,:][:,None]))\n",
    "    bottomedge_u = usol[-1,:][:,None]\n",
    "    \n",
    "    all_X_u_train = np.vstack([leftedge_x, rightedge_x, bottomedge_x, topedge_x])\n",
    "    all_u_train = np.vstack([leftedge_u, rightedge_u, bottomedge_u, topedge_u])  \n",
    "     \n",
    "    #choose random N_u points for training\n",
    "    idx = np.random.choice(all_X_u_train.shape[0], N_u, replace=False) \n",
    "    \n",
    "    X_u_train = all_X_u_train[idx[0:N_u], :] #choose indices from  set 'idx' (x,t)\n",
    "    u_train = all_u_train[idx[0:N_u],:]      #choose corresponding u\n",
    "    \n",
    "    '''Collocation Points'''\n",
    "\n",
    "    # Latin Hypercube sampling for collocation points \n",
    "    # N_f sets of tuples(x,t)\n",
    "    X_f = lb + (ub-lb)*lhs(2,N_f) \n",
    "    X_f_train = np.vstack((X_f, X_u_train)) # append training points to collocation points \n",
    "    \n",
    "    return X_f_train, X_u_train, u_train \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PINN \n",
    "\n",
    "$W \\in \\mathcal{R}^{n_{l-1}\\times{n_l}}$ \n",
    "\n",
    "Creating sequential layers using the $\\textit{class}$ tf.Module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Sequentialmodel(tf.Module): \n",
    "    def __init__(self, layers, name=None):\n",
    "\n",
    "        self.W = []  #Weights and biases\n",
    "        self.parameters = 0 #total number of parameters\n",
    "\n",
    "        for i in range(len(layers)-1):\n",
    "\n",
    "            input_dim = layers[i]\n",
    "            output_dim = layers[i+1]\n",
    "\n",
    "            #Xavier standard deviation \n",
    "            std_dv = np.sqrt((2.0/(input_dim + output_dim)))\n",
    "\n",
    "            #weights = normal distribution * Xavier standard deviation + 0\n",
    "            w = tf.random.normal([input_dim, output_dim], dtype = 'float64') * std_dv\n",
    "\n",
    "            w = tf.Variable(w, trainable=True, name = 'w' + str(i+1))\n",
    "\n",
    "            b = tf.Variable(tf.cast(tf.zeros([output_dim]), dtype = 'float64'), trainable = True, name = 'b' + str(i+1))\n",
    "\n",
    "            self.W.append(w)\n",
    "            self.W.append(b)\n",
    "\n",
    "            self.parameters +=  input_dim * output_dim + output_dim\n",
    "            \n",
    "        self.X = np.zeros(self.parameters) #store iterates\n",
    "        self.G = np.zeros(self.parameters) #store gradients\n",
    "        self.store = np.zeros((max_iter,2)) #store computed values for plotting\n",
    "        self.iter_counter = 0 # iteration counter for optimizer\n",
    "    \n",
    "    def evaluate(self,x):\n",
    "        \n",
    "        #preprocessing input \n",
    "        x = (x - lb)/(ub - lb) #feature scaling\n",
    "        \n",
    "        a = x\n",
    "        \n",
    "        for i in range(len(layers)-2):\n",
    "            \n",
    "            z = tf.add(tf.matmul(a, self.W[2*i]), self.W[2*i+1])\n",
    "            a = tf.nn.tanh(z)\n",
    "            \n",
    "        a = tf.add(tf.matmul(a, self.W[-2]), self.W[-1]) # For regression, no activation to last layer\n",
    "        \n",
    "        return a\n",
    "    \n",
    "    def get_weights(self):\n",
    "\n",
    "        parameters_1d = []  # [.... W_i,b_i.....  ] 1d array\n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "            \n",
    "            w_1d = tf.reshape(self.W[2*i],[-1])   #flatten weights \n",
    "            b_1d = tf.reshape(self.W[2*i+1],[-1]) #flatten biases\n",
    "            \n",
    "            parameters_1d = tf.concat([parameters_1d, w_1d], 0) #concat weights \n",
    "            parameters_1d = tf.concat([parameters_1d, b_1d], 0) #concat biases\n",
    "        \n",
    "        return parameters_1d\n",
    "        \n",
    "    def set_weights(self,parameters):\n",
    "                \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            shape_w = tf.shape(self.W[2*i]).numpy() # shape of the weight tensor\n",
    "            size_w = tf.size(self.W[2*i]).numpy() #size of the weight tensor \n",
    "            \n",
    "            shape_b = tf.shape(self.W[2*i+1]).numpy() # shape of the bias tensor\n",
    "            size_b = tf.size(self.W[2*i+1]).numpy() #size of the bias tensor \n",
    "                        \n",
    "            pick_w = parameters[0:size_w] #pick the weights \n",
    "            self.W[2*i].assign(tf.reshape(pick_w,shape_w)) # assign  \n",
    "            parameters = np.delete(parameters,np.arange(size_w),0) #delete \n",
    "            \n",
    "            pick_b = parameters[0:size_b] #pick the biases \n",
    "            self.W[2*i+1].assign(tf.reshape(pick_b,shape_b)) # assign \n",
    "            parameters = np.delete(parameters,np.arange(size_b),0) #delete \n",
    "\n",
    "            \n",
    "    def loss_BC(self,x,y):\n",
    "\n",
    "        loss_u = tf.reduce_mean(tf.square(y-self.evaluate(x)))\n",
    "        return loss_u\n",
    "\n",
    "    def loss_PDE(self, x_to_train_f):\n",
    "    \n",
    "        g = tf.Variable(x_to_train_f, dtype = 'float64', trainable = False)\n",
    "\n",
    "        k = 1    \n",
    "\n",
    "        x_1_f = g[:,0:1]\n",
    "        x_2_f = g[:,1:2]\n",
    "\n",
    "        with tf.GradientTape(persistent=True) as tape:\n",
    "\n",
    "            tape.watch(x_1_f)\n",
    "            tape.watch(x_2_f)\n",
    "\n",
    "            g = tf.stack([x_1_f[:,0], x_2_f[:,0]], axis=1)\n",
    "\n",
    "            u = self.evaluate(g)\n",
    "            u_x_1 = tape.gradient(u,x_1_f)\n",
    "            u_x_2 = tape.gradient(u,x_2_f)\n",
    "\n",
    "        u_xx_1 = tape.gradient(u_x_1,x_1_f)\n",
    "        u_xx_2 = tape.gradient(u_x_2,x_2_f)\n",
    "\n",
    "        del tape\n",
    "\n",
    "        q = -( (a_1*np.pi)**2 + (a_2*np.pi)**2 - k**2 ) * np.sin(a_1*np.pi*x_1_f) * np.sin(a_2*np.pi*x_2_f)\n",
    "\n",
    "        f = u_xx_1 + u_xx_2 + k**2 * u - q #residual\n",
    "\n",
    "        loss_f = tf.reduce_mean(tf.square(f))\n",
    "\n",
    "        return loss_f, f\n",
    "    \n",
    "    def loss(self,x,y,g):\n",
    "\n",
    "        loss_u = self.loss_BC(x,y)\n",
    "        loss_f, f = self.loss_PDE(g)\n",
    "\n",
    "        loss = loss_u + loss_f\n",
    "\n",
    "        return loss, loss_u, loss_f \n",
    "    \n",
    "    def optimizerfunc(self,parameters):\n",
    "        \n",
    "        self.set_weights(parameters)\n",
    "       \n",
    "        with tf.GradientTape() as tape:\n",
    "            tape.watch(self.trainable_variables)\n",
    "            \n",
    "            loss_val, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "            \n",
    "        grads = tape.gradient(loss_val,self.trainable_variables)\n",
    "                \n",
    "        del tape\n",
    "        \n",
    "        grads_1d = [ ] #store 1d grads \n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            grads_w_1d = tf.reshape(grads[2*i],[-1]) #flatten weights \n",
    "            grads_b_1d = tf.reshape(grads[2*i+1],[-1]) #flatten biases\n",
    "\n",
    "            grads_1d = tf.concat([grads_1d, grads_w_1d], 0) #concat grad_weights \n",
    "            grads_1d = tf.concat([grads_1d, grads_b_1d], 0) #concat grad_biases\n",
    "        \n",
    "        return loss_val.numpy(), grads_1d.numpy()\n",
    "    \n",
    "    def optimizer_callback(self,parameters):\n",
    "                \n",
    "        loss_value, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "        \n",
    "        u_pred = self.evaluate(X_u_test)\n",
    "        error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)\n",
    "        \n",
    "        alpha, count, f_newval = self.LbfgsInvHessProduct(parameters)\n",
    "        \n",
    "        # if line search fails, it returns None, but difficult to handle while plotting\n",
    "        # so change it to 0\n",
    "        \n",
    "        if alpha == None:\n",
    "            alpha = 0\n",
    "        \n",
    "        tf.print(loss_value, f_newval, (loss_value-f_newval)/loss_value, alpha, count)\n",
    "        \n",
    "    def LbfgsInvHessProduct(self,parameters):\n",
    "\n",
    "        self.iter_counter += 1  #update iteration counter \n",
    "\n",
    "        x_k = parameters  \n",
    "\n",
    "        self.X = np.vstack((x_k.T,self.X)) #stack latest value on top row\n",
    "\n",
    "        old_fval,g_k = self.optimizerfunc(parameters) #obtain grads and loss value\n",
    "\n",
    "        self.G = np.vstack((g_k.T,self.G)) #stack latest grads on top row\n",
    "\n",
    "        n_corrs = min(self.iter_counter, maxcor) #for iterations < maxcor, we will take all available updates\n",
    "        \n",
    "        sk = self.X = self.X[:n_corrs] #select top 'n_corrs' x values, with latest value on top by construction\n",
    "        yk = self.G = self.G[:n_corrs] #select top 'n_corrs' gradient values, with latest value on top by construction \n",
    "\n",
    "        #linear operator B_k_inv    \n",
    "        hess_inv = scipy.optimize.LbfgsInvHessProduct(sk,yk) #instantiate class\n",
    "\n",
    "        p_k = - hess_inv.matvec(g_k) #p_k = -B_k_inv * g_k\n",
    "\n",
    "        gkpk = np.dot(p_k,g_k) #term 1 in report\n",
    "\n",
    "        norm_p_k_sq = (np.linalg.norm(p_k,ord=2))**2 # norm squared\n",
    "               \n",
    "        #store the values\n",
    "        self.store[self.iter_counter-1] = [gkpk,norm_p_k_sq]\n",
    "        \n",
    "        def line_search_armijo(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1):\n",
    "            \n",
    "            \"\"\"\n",
    "            Minimize over alpha, the function ``f(xk+alpha pk)``\n",
    "            \n",
    "            Parameters\n",
    "            ----------\n",
    "            f : callable\n",
    "                Function to be minimized.\n",
    "            xk : array_like\n",
    "                Current point.\n",
    "            pk : array_like\n",
    "                Search direction.\n",
    "            gfk : array_like\n",
    "                Gradient of `f` at point `xk`.\n",
    "            old_fval : float\n",
    "                Value of `f` at point `xk`.\n",
    "            args : tuple, optional\n",
    "                Optional arguments.\n",
    "            c1 : float, optional\n",
    "                Value to control stopping criterion.\n",
    "            alpha0 : scalar, optional\n",
    "                Value of `alpha` at start of the optimization.\n",
    "                \n",
    "            Returns\n",
    "            -------\n",
    "            alpha\n",
    "            f_count\n",
    "            f_val_at_alpha\n",
    "            Notes\n",
    "            -----\n",
    "            Uses the interpolation algorithm (Armijo backtracking) as suggested by\n",
    "            Wright and Nocedal in 'Numerical Optimization', 1999, pp. 56-57\n",
    "            \"\"\"\n",
    "            \n",
    "            xk = np.atleast_1d(xk)\n",
    "            fc = [0]\n",
    "\n",
    "            def phi(alpha1):\n",
    "                fc[0] += 1\n",
    "                return f(xk + alpha1*pk, *args)\n",
    "\n",
    "            if old_fval is None:\n",
    "                phi0 = phi(0.)\n",
    "            else:\n",
    "                phi0 = old_fval  # compute f(xk) -- done in past loop\n",
    "\n",
    "            derphi0 = np.dot(gfk, pk)\n",
    "            if derphi0 > 0:\n",
    "                print('Warning dephi0 :' + str(derphi0))\n",
    "            alpha, phi1 = scalar_search_armijo(phi, phi0, derphi0, c1=c1,\n",
    "                                               alpha0=alpha0)\n",
    "            return alpha, fc[0], phi1\n",
    "\n",
    "\n",
    "        def line_search_BFGS(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1):\n",
    "            \n",
    "            \"\"\"\n",
    "            Compatibility wrapper for `line_search_armijo`\n",
    "            \"\"\"\n",
    "            \n",
    "            r = line_search_armijo(f, xk, pk, gfk, old_fval, args=args, c1=c1,\n",
    "                                   alpha0=alpha0)\n",
    "            return r[0], r[1], 0, r[2]\n",
    "\n",
    "\n",
    "        def scalar_search_armijo(phi, phi0, derphi0, c1=1e-4, alpha0=1, amin=0):\n",
    "            \n",
    "            \"\"\"\n",
    "            Minimize over alpha, the function ``phi(alpha)``.\n",
    "            Uses the interpolation algorithm (Armijo backtracking) as suggested by\n",
    "            Wright and Nocedal in 'Numerical Optimization', 1999, pp. 56-57\n",
    "            alpha > 0 is assumed to be a descent direction.\n",
    "            \n",
    "            Returns\n",
    "            -------\n",
    "            alpha\n",
    "            phi1\n",
    "            \n",
    "            \"\"\"\n",
    "            \n",
    "            phi_a0 = phi(alpha0)\n",
    "            if (phi_a0 <=  phi0 + c1*alpha0*derphi0):\n",
    "                return alpha0, phi_a0\n",
    "\n",
    "            # Otherwise, compute the minimizer of a quadratic interpolant:\n",
    "\n",
    "            alpha1 = -(derphi0) * alpha0**2 / 2.0 / (phi_a0 - phi0 - derphi0 * alpha0)\n",
    "            phi_a1 = phi(alpha1)\n",
    "\n",
    "            if (phi_a1 <= phi0 + c1*alpha1*derphi0):\n",
    "                return alpha1, phi_a1\n",
    "\n",
    "            # Otherwise, loop with cubic interpolation until we find an alpha which\n",
    "            # satisfies the first Wolfe condition (since we are backtracking, we will\n",
    "            # assume that the value of alpha is not too small and satisfies the second\n",
    "            # condition.\n",
    "\n",
    "            while alpha1 > amin:       # we are assuming alpha>0 is a descent direction\n",
    "                factor = alpha0**2 * alpha1**2 * (alpha1-alpha0)\n",
    "                a = alpha0**2 * (phi_a1 - phi0 - derphi0*alpha1) - \\\n",
    "                    alpha1**2 * (phi_a0 - phi0 - derphi0*alpha0)\n",
    "                a = a / factor\n",
    "                b = -alpha0**3 * (phi_a1 - phi0 - derphi0*alpha1) + \\\n",
    "                    alpha1**3 * (phi_a0 - phi0 - derphi0*alpha0)\n",
    "                b = b / factor\n",
    "\n",
    "                alpha2 = (-b + np.sqrt(abs(b**2 - 3 * a * derphi0))) / (3.0*a)\n",
    "                phi_a2 = phi(alpha2)\n",
    "\n",
    "                if (phi_a2 <= phi0 + c1*alpha2*derphi0):\n",
    "                    return alpha2, phi_a2\n",
    "\n",
    "                if (alpha1 - alpha2) > alpha1 / 2.0 or (1 - alpha2/alpha1) < 0.96:\n",
    "                    alpha2 = alpha1 / 2.0\n",
    "\n",
    "                alpha0 = alpha1\n",
    "                alpha1 = alpha2\n",
    "                phi_a0 = phi_a1\n",
    "                phi_a1 = phi_a2\n",
    "\n",
    "            # Failed to find a suitable step length\n",
    "            return None, phi_a1\n",
    "        \n",
    "        # wrapper for line_search_BFGS\n",
    "        def ls_function(x):\n",
    "            val, _ = self.optimizerfunc(x)\n",
    "            return val\n",
    "        \n",
    "        alpha, count, _, f_newval = line_search_BFGS(ls_function, x_k, p_k, g_k, old_fval, args=(), c1=1e-4, alpha0=1)    \n",
    "        \n",
    "        return alpha, count, f_newval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Loss Function*\n",
    "\n",
    "The loss function consists of two parts:\n",
    "1. **loss_BC**: MSE error of boundary losses\n",
    "2. **loss_PDE**: MSE error of collocation points satisfying the PDE\n",
    "\n",
    "**loss** = loss_BC + loss_PDE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_u = 400 #Total number of data points for 'u'\n",
    "N_f = 10000 #Total number of collocation points \n",
    "\n",
    "# Training data\n",
    "X_f_train, X_u_train, u_train = trainingdata(N_u,N_f)\n",
    "\n",
    "layers = np.array([2, 50, 50, 50, 1]) #3 hidden layers\n",
    "\n",
    "maxcor = 200 \n",
    "max_iter = 5000\n",
    "\n",
    "PINN = Sequentialmodel(layers)\n",
    "\n",
    "init_params = PINN.get_weights().numpy()\n",
    "\n",
    "start_time = time.time() \n",
    "\n",
    "# train the model with Scipy L-BFGS optimizer\n",
    "results = scipy.optimize.minimize(fun = PINN.optimizerfunc, \n",
    "                                  x0 = init_params, \n",
    "                                  args=(), \n",
    "                                  method='L-BFGS-B', \n",
    "                                  jac= True,        # If jac is True, fun is assumed to return the gradient along with the objective function\n",
    "                                  callback = PINN.optimizer_callback, \n",
    "                                  options = {'disp': None,\n",
    "                                            'maxcor': maxcor, \n",
    "                                            'ftol': 1 * np.finfo(float).eps,  #The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol\n",
    "                                            'gtol': 5e-10, \n",
    "                                            'maxfun':  max_iter*1.5, \n",
    "                                            'maxiter': max_iter,\n",
    "                                            'iprint': -1,   #print update every 50 iterations\n",
    "                                            'maxls': 50})\n",
    "\n",
    "elapsed = time.time() - start_time                \n",
    "print('Training time: %.2f' % (elapsed))\n",
    "\n",
    "print(results)\n",
    "\n",
    "# np.savetxt('values_stiff.txt', PINN.store)\n",
    "\n",
    "PINN.set_weights(results.x)\n",
    "\n",
    "''' Model Accuracy ''' \n",
    "u_pred = PINN.evaluate(X_u_test)\n",
    "\n",
    "error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)        # Relative L2 Norm of the error (Vector)\n",
    "print('Test Error: %.5f'  % (error_vec))\n",
    "\n",
    "u_pred = np.reshape(u_pred,(256,256),order='F') \n",
    "\n",
    "# #Residual plot\n",
    "# loss_f, f_plot = PINN.loss_PDE(X_u_test)\n",
    "# plt.scatter(X_u_test[:,0:1], X_u_test[:,1:2], c=f_plot, cmap = 'jet')\n",
    "# plt.axis('scaled')\n",
    "# plt.colorbar()\n",
    "# plt.savefig('Stiff_Helmholtz_residual.png', dpi = 500)\n",
    "\n",
    "# #plot gkpk\n",
    "# plt.semilogy(PINN.store[:,0])\n",
    "# plt.yscale('symlog')\n",
    "# plt.savefig('gkpk_stiff.png', dpi = 500)\n",
    "\n",
    "# #plot norm_p_k_sq\n",
    "# plt.semilogy(PINN.store[:,1])\n",
    "# plt.yscale('symlog')\n",
    "# plt.savefig('norm_p_k_sq_stiff.png', dpi = 500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = np.loadtxt('prove_stiffness/non_stiff.txt', comments = \"#\", delimiter = \" \", unpack = False)\n",
    "s2 = np.loadtxt('prove_stiffness/stiff.txt', comments = \"#\", delimiter = \" \", unpack = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison of actual loss and predicted loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Non-stiff case prediction error:0.0038351804643316583\n",
      "Stiff case prediction error:4.1792986909953895e-05\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABC8klEQVR4nO3deVzVVf748dfhsoOCyCKLCLkloqKiZi5puVVqaZa2TabWOGU5NpXV9G37TZNTzVQ6lVqWs5TamKaWmpoaalpuuG+oqCgIsik73Ht+f3wuiMqqXC7L+/l43Af3nvu5n8/7IPLmfM6mtNYIIYQQ5XGwdwBCCCHqNkkUQgghKiSJQgghRIUkUQghhKiQJAohhBAVcrR3ALbg6+urw8LC7B2GEELUKzt37rygtfa7urxBJoqwsDB27Nhh7zCEEKJeUUqdKqtcbj0JIYSoUL1IFEopD6XUDqXUcHvHIoQQjY1dEoVS6gulVLJSav9V5cOUUkeUUnFKqZdKvTUd+KZ2oxRCCAH266OYD/wT+HdxgVLKBHwMDAYSgO1KqeVAMHAQcK39MIVonAoLC0lISCAvL8/eoQgbcHV1JSQkBCcnpyodb5dEobWOUUqFXVXcE4jTWp8AUEotBO4BPAEPIALIVUqt1Fpbrj6nUupJ4EmA0NBQG0YvRMOXkJBAkyZNCAsLQyll73BEDdJak5qaSkJCAuHh4VX6TF0a9RQMnCn1OgHopbWeAqCUGg9cKCtJAGit5wJzAaKjo2WlQyFuQF5eniSJBkopRfPmzUlJSanyZ+pSoqiQ1nq+vWMQojGRJNFwVfffti6NejoLtCz1OsRaVmsy8jL4+y9/5+t9XyPLrwshhKEutSi2A22VUuEYCWIc8FBtBpCRl8GMLTO4kHOB3Ym7eW/Ie7V5eSGEqJPsNTx2AbAVaK+USlBKTdRaFwFTgB+BQ8A3WusDtRlXmHcY8VPjGRg2kH9u/ydnMs9U/iEhhM189913KKU4fPhwpcd++OGH5OTkXPe15s+fz5QpU8qN46233gLgjTfewN3dneTk5JL3PT09r/u6VXF13e666y4yMjIAmDlzJh06dODhhx8mPz+fQYMGERUVxaJFixg3bhzHjh274evbJVForR/UWgdqrZ201iFa63nW8pVa63Za69Za67ftEZuHswcfDP2AQnMhf9n0F3uEIISwWrBgAX379mXBggWVHnujiaIi7777Lk899VTJa19fX/7+97/b5FplubpuK1euxNvbG4BPPvmEtWvX8tVXX7F7924AYmNjGTt2LH/4wx949913b/j6damPos7o0qILd7a9kxVHVlBgLrB3OELYV/YZuHi0Zh/ZlbfWs7Ky2Lx5M/PmzWPhwoUl5Wazmeeff57IyEg6d+7MrFmzmDlzJufOnWPgwIEMHDgQuPKv/MWLFzN+/HgAVqxYQa9evejatSuDBg3i/PnzFcZx9OhRXFxc8PX1LSmbMGECixYtIi0t7Zrj//GPfxAZGUlkZCQffvghAPHx8XTo0IEnnniCjh07MmTIEHJzc6/9Vmdnc/fdd9OlSxciIyNZtGhRmXULCwvjwoULTJ48mRMnTnDnnXfyt7/9jUceeYTt27cTFRXF8ePH6devH+vWraOoqKjS73dFJFGU453b3+GzEZ+RlnvtD4IQwvaWLVvGsGHDaNeuHc2bN2fnzp0AzJ07l/j4eGJjY9m7dy8PP/wwzz77LEFBQWzYsIENGzZUeN6+ffuybds2du/ezbhx4yr9i3vLli1069btijJPT08mTJjARx99dEX5zp07+fLLL/n111/Ztm0bn332Wclf+ceOHePpp5/mwIEDeHt78+23315zrdWrVxMUFMSePXvYv38/w4YNq7Bus2fPLnlv+vTpfP755/Tr14/Y2Fhat26Ng4MDbdq0Yc+ePRXWsTJ1qTO7TokMiMTJ5MSFnAu08Gxh73CEsB+PlpUfYwMLFixg6tSpAIwbN44FCxbQvXt31q1bx+TJk3F0NH59+fj4VOu8CQkJjB07lsTERAoKCiqddJaYmIif3zUrb/Pss88SFRXF888/X1K2efNmRo0ahYeHBwCjR49m06ZNjBw5kvDwcKKiogDo3r078fHx15yzU6dO/OlPf2L69OkMHz6cfv36VatuZfH39+fcuXN07979us8hLYoKpOel88iSR9h4cqO9QxGiUUlLS2P9+vVMmjSJsLAw3nvvPb755ptqDVsvPVeg9FIkzzzzDFOmTGHfvn3MmTOn0mVK3NzcyjzG29ubhx56iI8//rhK8bi4uJQ8N5lMFBUVcebMGaKiooiKimL27Nm0a9eOXbt20alTJ1599dWSDvQbkZeXh5ub2w2dQxJFBdr5tOPwhcN8uedLe4ciRKOyePFiHn30UU6dOkV8fDxnzpwhPDycTZs2MXjwYObMmVNy3724n6BJkyZcunSp5BwBAQEcOnQIi8XC0qVLS8ozMzMJDg4G4F//+lelsXTo0IG4uLgy33vuueeuiKVfv35899135OTkkJ2dzdKlSytsFbRs2ZLY2FhiY2OZPHky586dw93dnUceeYQXXniBXbt2lVm36jh69CiRkZHX9dlikigq4OPuw+3ht/PD0R/IK5TF0YSoLQsWLGDUqFFXlN13330sWLCASZMmERoaSufOnenSpQtff/01AE8++STDhg0r6fCdMWMGw4cP59ZbbyUwMLDkPG+88Qb3338/3bt3v6KDujz9+/dn9+7dZbZmfH19GTVqFPn5+QB069aN8ePH07NnT3r16sWkSZPo2rVrleu9b98+evbsSVRUFG+++SavvvpqmXWrqvPnz+Pm5kaLFjd2+1w1xBnI0dHRuqZ2uFt2ZBn3LryXJ7o9wdwRc2vknELUdYcOHaJDhw72DqPOmDp1KiNGjGDQoEH2DqVaPvjgA5o2bcrEiROvea+sf2Ol1E6tdfTVx0qLohL3tL+H33X5HZ/t+ox/x/678g8IIRqcV155xWZzNGzJ29ubxx577IbPI6OequDzEZ/Tr2U/OrfojEVbcFCSX4VoTAICAhg5cqS9w6i2xx9/vEbOI7/xqsDJ5MS4TuMwW8ycz6p4co4QQjQ0kiiqyNPZk5hTMfSe15vk7OTKPyCEEA2EJIpq6Bnck7OXzjJp+SR7hyKEELVGEkU19Antw5QeU1hxdAVrT6y1dzhCNGgmk4moqCgiIyO5//77b6gzefz48SxevBiASZMmcfDgwXKP3bhxI7/88ku1r1G8/lJDJImiml6/7XX83P2Yvna6bG4khA25ubkRGxvL/v37cXZ2Zvbs2Ve8f70L3X3++edERESU+/71JoqGTBJFNXm7efNinxfZnbSbDfEVLz4mhKgZ/fr1Iy4ujo0bN9KvXz9GjhxJREQEZrOZF154gR49etC5c2fmzJkDgNaaKVOm0L59ewYNGnTF3hEDBgygeJ7V6tWr6datG126dOGOO+4gPj6e2bNn88EHHxAVFcWmTZtISUnhvvvuo0ePHvTo0YMtW7YAkJqaypAhQ+jYsSOTJk1q0H84yvDY6zC111TCvMLw9/BHay17C4uGb92Aa8tCH4B2T0FRDmy869r3bxpvPPIuwOYxV743aGOVL11UVMSqVasYNmwYALt27WL//v2Eh4czd+5cvLy82L59O/n5+fTp04chQ4awe/dujhw5wsGDBzl//jwRERFMmDDhivOmpKTwxBNPEBMTQ3h4OGlpafj4+DB58mQ8PT1LFvt76KGHmDZtGn379uX06dMMHTqUQ4cO8eabb9K3b19ee+01fvjhB+bNm1flOtU3kiiug5PJicGtBxOXFseFnAv4eVy7sqQQ4sbk5uaWrLbar18/Jk6cyC+//ELPnj1LVnxds2YNe/fuLel/yMzM5NixY8TExPDggw9iMpkICgri9ttvv+b827Zto3///iXnKm8V2nXr1l3Rp3Hx4kWysrKIiYlhyZIlANx99900a9asxupe10iiuE5erl7M2DKDsxfPsmXCFmlViIatohaAo3vF77v6VqsFUay4j+JqxUt4g3GLadasWQwdOvSKY1auXFnt65XHYrGwbds2XF1da+yc9Y30UdyAHkE92JqwlU+2f2LvUIRolIYOHcqnn35KYWEhYKyUmp2dTf/+/Vm0aBFms5nExMQyNzO65ZZbiImJ4eTJk0D5q9AOGTKEWbNmlbwuTl79+/cvWZBw1apVpKen26SOdYEkihvw/K3P07VFV/68/s+cSD9h73CEaHQmTZpEREQE3bp1IzIykt///vcUFRUxatQo2rZtS0REBL/73e/o3bv3NZ/18/Nj7ty5jB49mi5dujB27FgARowYwdKlS0s6s2fOnMmOHTvo3LkzERERJaOvXn/9dWJiYujYsSNLliwhNDS0Vutem2T12Bu0O3E3fb/sSxufNmyftB1nR+daua4QtiSrxzZ8snpsLeoa2JVZd87iRPoJGS4rhGiQpDO7BkzoOoHugd0pshSRlpuGj1v19vAVQoi6TFoUNaRzQGfcnNz466a/ciz1mL3DEUKIGiOJooYopXAxufDJ9k/4vw3/Z+9whBCixkiiqEGtfVoztPVQvj/6PReyG+biYEKIxkcSRQ17qd9LZBdm87ctf7N3KEIIUSMkUdSwXsG9uCP8DubumiutCiFEgyCJwgbevv1twrzCOJhS/pr3QojKfffddyilOHz4cKXHfvjhhze0Z8X8+fOZMmVKuXG89dZbABw5coQBAwYQFRVFhw4dePLJJwFjxnbppUOWL1/OjBkzAGMBwl69etG1a1c2bdrE//73Pzp06MDAgQPZt28f48ePv+64a4MkChvoFdKLFQ+uwMPZg0v5lyr/gBCiTAsWLKBv374sWLCg0mNvNFFU5N133+Wpp54C4Nlnn2XatGnExsZy6NAhnnnmGeDaRDFy5EheeuklAH766Sc6derE7t276devH/PmzeOzzz5jw4YNdOrUiYSEBE6fPm2T2GuCzKOwkeCmwSRcSmDC8gn8Y8g/aOnV0t4hCXFdzmSeIbcot0bP6eboVun/iaysLDZv3syGDRsYMWIEb775JgBms5np06ezevVqHBwceOKJJ9Bac+7cOQYOHIivry8bNmzA09OTrKwsABYvXsz333/P/PnzWbFiBX/5y18oKCigefPmfPXVVwQEBJQbx9GjR3FxccHX1xeAxMREQkJCSt7v1KkTBQUFvPbaa+Tm5rJ582ZefvllcnNz2bFjB5MmTeLFF18seT1q1Cg2b97MxIkTGTlyJO+99x4jRoxg4cKFvPjiizf6rbUJaVHYiMnBhIvJhe+PfM/4ZeMb9KYmQtjCsmXLGDZsGO3ataN58+bs3LkTgLlz5xIfH09sbCx79+7l4Ycf5tlnnyUoKIgNGzaUuQBgaX379mXbtm3s3r2bcePG8e6771Z4/JYtW+jWrVvJ62nTpnH77bdz55138sEHH5CRkYGzszNvvfUWY8eOJTY2tmTdKICoqKgr3nv99deJjo7mq6++4r333gMgOjqaTZs2Xe+3yuakRWFD3YO6M73vdN78+U3mx87n8a6P2zskIarNXq3hBQsWMHXqVADGjRvHggUL6N69O+vWrWPy5Mk4Ohq/vsrbR6I8CQkJjB07lsTERAoKCkr2oyhPYmIifn6X95x5/PHHGTp0KKtXr2bZsmXMmTOHPXv2VLN2V/L39+fcuXM3dA5bqvMtCqXUvUqpz5RSi5RSQ+wdT3W92v9VIvwieHHdi6Rkp9g7HCHqhbS0NNavX8+kSZMICwvjvffe45tvvqlWy7z0HjF5eXklz5955hmmTJnCvn37mDNnzhXvlcXNze2aY4KCgpgwYQLLli3D0dGR/fv3VzmusuTl5eHm5nZD57AluyQKpdQXSqlkpdT+q8qHKaWOKKXilFIvAWitv9NaPwFMBsaWdb66zNHBkc9GfEZ6bjovrqub9x+FqGsWL17Mo48+yqlTp4iPj+fMmTOEh4ezadMmBg8ezJw5cygqKgLK30ciICCAQ4cOYbFYWLp0aUl5ZmYmwcHBAPzrX/+qNJYOHToQFxdX8nr16tUl+18kJSWRmppKcHDwNdevjqNHjxIZGXldn60N9mpRzAeGlS5QSpmAj4E7gQjgQaVURKlDXrW+X+/c2vJW3h/yPuMix5FflG/vcISo8xYsWMCoUaOuKLvvvvtYsGABkyZNIjQ0lM6dO9OlS5eSzYOefPJJhg0bxsCBAwGYMWMGw4cP59ZbbyUwMLDkPG+88Qb3338/3bt3L+mgrkj//v3ZvXt3SWtmzZo1REZG0qVLF4YOHcp7771HixYtGDhwIAcPHiQqKopFixZVq74bNmzg7rvvrtZnapPd9qNQSoUB32utI62vewNvaK2HWl+/bD10hvWxVmu9roLzPQk8CRAaGtr91KlTNoy++grMBRxMOYijgyPtfNrJvhWiTpP9KK40depURowYwaBBg2r83Pn5+dx2221s3ry5pN+lNtTX/SiCgTOlXidYy54BBgFjlFKTy/uw1nqu1jpaax1duuOpWvJSYNcLcLHmV391NjkT6BnIY989JreghKhnXnnlFZvN0Th9+jQzZsyo1SRRXXUpUZRJaz1Ta91daz1Zaz3bphezFMKxj2HfGzY5fYBnAC2btmT2jtkcTJZZ20LUFwEBAYwcOdIm527bti0DBgywyblrSl1KFGeB0uPwQqxltcc9CPz6wfn1YKNbcv8Y8g9cHF2YuGKizK0QQtQLdSlRbAfaKqXClVLOwDhgea1HETgU8pIg48aGu5XnJp+beLHPi2xL2Ma3B7+1yTWEEKIm2Wt47AJgK9BeKZWglJqotS4CpgA/AoeAb7TWB2o9OB9rP07mXptd4vnezxPqFcqcnXOwaIvNriOEEDXBLolCa/2g1jpQa+2ktQ7RWs+zlq/UWrfTWrfWWr9tj9jw7mx8tUGHdjEXRxf+M+o/vH7b6yRlJdnsOkLUZyaTiaioKCIjI7n//vtvqDN5/PjxLF68GIBJkyZx8GD5fYQbN27kl19+qfY1wsLCuHChYW4tUJduPdUNLt4waAuEP2LTy/Rv1Z/gpsGcyTzD+azzNr2WEPWRm5sbsbGx7N+/H2dnZ2bPvnIsS/GEu+r6/PPPiYiIKPf9600UDZkkirK4B0GRbYbClebn4ceY/41h2o/TbH4tIeqzfv36ERcXx8aNG+nXrx8jR44kIiICs9nMCy+8QI8ePejcuTNz5swBQGvNlClTaN++PYMGDSI5ObnkXAMGDGDHjh2AMcu6W7dudOnShTvuuIP4+Hhmz57NBx98QFRUFJs2bSIlJYX77ruPHj160KNHD7Zs2QJAamoqQ4YMoWPHjkyaNKlBD06puwN37enCVoj/D/RfDg62+xZ5Onsy+KbBzI+dz9iOY7nn5ntsdi0hbsSA+QOuKXug4wM81eMpcgpzuOuru655f3zUeMZHjedCzgXGfDPmivc2jt9Y5WsXFRWxatUqhg0zFnPYtWsX+/fvJzw8nLlz5+Ll5cX27dvJz8+nT58+DBkyhN27d3PkyBEOHjzI+fPniYiIYMKECVecNyUlhSeeeIKYmBjCw8NJS0vDx8eHyZMn4+npyfPPPw/AQw89xLRp0+jbty+nT59m6NChHDp0iDfffJO+ffvy2muv8cMPPzBv3rwq16m+kURRltxEOLcKLh0Hr/Y2vdRHwz7i51M/M2nFJKKDogluGmzT6wlRX+Tm5hIVFQUYLYqJEyfyyy+/0LNnz5IVX9esWcPevXtL+h8yMzM5duwYMTExPPjgg5hMJoKCgrj99tuvOf+2bdvo379/ybnKW4V23bp1V/RpXLx4kaysLGJiYliyZAkAd999N82aNauxutc1kijK4t3R+Jp5wOaJoolLExbet5A+X/Rh4vKJrH5ktU2vJ8T1qKgF4O7kXuH7vu6+1WpBFCvuo7iah4dHyXOtNbNmzWLo0KFXHFN6p7kbZbFY2LZtG66urjV2zvpG+ijK0sSaHLKO18rlegT34Jmez5CWm0ZGbkatXFOIhmDo0KF8+umnJau5Hj16lOzsbPr378+iRYswm80kJiaWuZnRLbfcQkxMDCdPngTKX4V2yJAhzJo1q+R1cfLq379/yYKEq1atIj093SZ1rAukRVEWj1Awudts0l1Z3hn0DgeSD5Cam4q3m3etXVeI+mzSpEnEx8fTrVs3tNb4+fnx3XffMWrUKNavX09ERAShoaH07t37ms/6+fkxd+5cRo8ejcViwd/fn7Vr1zJixAjGjBnDsmXLmDVrFjNnzuTpp5+mc+fOFBUV0b9/f2bPns3rr7/Ogw8+SMeOHbn11lsJDQ21w3egdtht9Vhbio6O1sWjGq7bugHg3Az6L6300JpyPus8O87tIDM/k4c6PVRr1xXiarJ6bMNXndVjpUVRnt7/hfxk0BZQtXOHzt/Dn5m/zSTmVAxdWnSho1/HWrmuEEJURPooyuPUBNBQlF1rl1RK8fFdH+NscmbE1yNk61QhRJ0giaI8Dq7w2+9h//+r1cu28WnD16O/Jj4jnud+fK5Wry1EaQ3xtrQwVPffVhJFeRxdwJwHyZtq/dJ3t7ubCV0n8N99/+XrvV/X+vWFcHV1JTU1VZJFA6S1JjU1tVrDfaWPoiLeneHcD8beFErV6qVn3TkLHzcfwpuFU2QpwtGGM8SFuFpISAgJCQmkpMjtz4bI1dWVkJCQKh8vv30q0iwKTi2AS3HQtG2tXtrNyY23Br7FgeQDnM44zU0+N9Xq9UXj5uTkVDJjWQi59VSR5j2Mr2nb7XJ5V0dXsgqyGPCvAfxnz3/sEoMQQkiiqIhPNATcAcp+Da9ugd1o5tqMp1c+TWxSrN3iEEI0XpIoKuLUBLq9D7632C2EJi5N+Pq+r3FQDkxYNkE6F4UQtU4SRWUcnI3RT3bU0b8jrw94nd1Juxm7eCxF5uvbsEUIIa6HdGZXZtdzkHsO7rLdHtpVMbXXVE6kneDQhUOcyjxFa5/Wdo1HCNF4SKKojGMTyLX/VqUOyoGZd87kZPpJ0vPSOZRyiJt9b0bV8rBdIUTjI7eeKuMeDAWpUAdu9yilCGsWhsnBxOD/DOahbx/Coi32DksI0cBJoqiMWzBoM+Qm2DsSwGhZRPpFMurmUSw8sJBHljwifRZCCJuSW0+V8WhpfM0+BZ5hdg2lmKPJkZl3zsSszXy641POXTrH9w99j6ezp71DE0I0QNKiqIxXJ2j1IDh6VH5sLVJK8cndn/DGbW/w86mfeXrl03IbSghhE9KiqEzTm6HDn8Ct6uui1KbXB7xOB78OBDcJ5siFI4R5h+Hm5GbvsIQQDYi0KCrjYAJzAeQl2zuScj3Q8QE6+nckMy+TWz6/hSWHltg7JCFEAyKJoio2jYb9b9g7igp5u3rT1KUpGfkZ3PfNfTz+3eNyK0oIUSMkUVSFawDkJto7ikpF+Eew5/d7GNdxHPP3zOfWebdyPO24vcMSQtRzkiiqwi0Q8uw/6a4qvN28+fq+r5lxxwx2J+3mhbUvkJabZu+whBD1mHRmV4VbEKRssssGRtdDKcX0vtMZ0W4EuUW5nEw/yd7zewnwCKCDXwd7hyeEqGckUVSFWxAUZUNBBrg0s3c0VRbhH4HWmvPZ57l13q0kZSXxdI+nefv2t3F2dLZ3eEKIekJuPVVFi0HQ/o9gKbB3JNWmlKKFZwu+f+h7OgV04v2t79Pyw5b887d/kl+Ub+/whBD1gCSKqmjeA8IfAZOLvSO5bhF+Efw66VcW3rcQf3d/nln1DF/u/pLMvEzZ40IIUaE6f+tJKeUBfAIUABu11l/VfhQmyD4NDi7g7F37l69BYyPHck/7e4g5FYO/pz9xaXG8s/kdmro05e3b3yawSaC9QxRC1DF2aVEopb5QSiUrpfZfVT5MKXVEKRWnlHrJWjwaWKy1fgIYWevBAuhCYy5FvB1ylA24OrkypM0QugR0obVPawothXwZ+yXhH4UzbfU0jqYetXeIQog6xF63nuYDw0oXKKVMwMfAnUAE8KBSKgIIAc5YDzPXYoyXOXmCo2e9mEtRHUopvF29WTZuGb9O+pWugV358NcP6flZT2KTYimyyKq0Qgg73XrSWscopcKuKu4JxGmtTwAopRYC9wAJGMkiFnv2qbgGGDvdNVA9g3uydeJWdp7byTOrniG/KJ995/ex7Mgy/N39GRg+UIbWCtFI1aU+imAutxzASBC9gJnAP5VSdwMryvuwUupJ4EmA0NDQmo/OtUW9mXR3I7oHdeeXib+QU5jD+azzfL3va46lHQOgV3Avnuj2BBO6TpCd9YRoROpSoiiT1jobeLwKx80F5gJER0fX/DAet0BIiavx09ZV7k7uhDcL58BTB9iTtIclh5cwb/c8Jq2YRGJWIhO7TqSFZwtJGEI0AnUpUZwFWpZ6HWItqxtuGg/Ne4KlCBzq0rfNtpxMTkQHRxMdHM1fbv8Lf//l77Rr3o5zl87x373/5WTGSUa0G0H/Vv3xcK5be3YIIWpGXfqNtx1oq5QKx0gQ44CH7BtSKX59wK2FMemuESWK0hyUAy/0eQGA1JxUjqQeYX7sfD7d8Sn+Hv58MfIL+rTsg7ebt30DFULUKHsNj10AbAXaK6USlFITtdZFwBTgR+AQ8I3W+oA94iuTOQ9St0Nukr0jqROauzfn85Gfk/x8Mp/e/Sk5hTkMXzCcexbdw9mLZ2XWtxANiGqIs3Kjo6P1jh07avakKVth7a1wy7/hpkdr9twNQFJWEutPrqeDbwcs2sKBlAN8e+hbxnYcy8h2I/F0kf28hajrlFI7tdbRV5c3znso18PdOpIq50zFxzVSLTxb8FAn405hgbmA2Ttns+nUJpYfWY6DcmBA2ADuanMXU2+ZimMjvXUnRH0l/2Oryj0IHD3g4mF7R1LnOZuc+WzEZ3w49EMWH1zM+vj1rDy2kryiPG5rdRvert609GqJi2P9XTtLiMakSonCut5SrtbaopRqB9wMrNJaF9o0urpEKWjSDi4esnck9YaHswePRT3GY1GPAZB4KRGzNrPy2EpeWvcSAZ4BfDD0A4a1GVbJmYQQ9lTVzuwYwFUpFQysAR7FWIajcfGKgItHwNx48mNNCmwSSEjTEEKbhnJ/xP1kF2Qz/OvhRH4SyR++/wOX8i/ZO0QhRBmqeutJaa1zlFITgU+01u8qpWJtGFfd1PYpCB4B5hwwedk7mnprWNthDGs7jHOXzjFj8wx+O/sb289t58iFIziZnPjx+I9E+EYwtM1QnExO9g5XiEavyolCKdUbeBiYaC0z2SakOqx5L2NPiqJscJZEcaOCmgQx886ZAOQW5pJVkMW5S+d45adXMGszQZ5BRAdH0zukN6M7jKaNTxsclGyhIkRtq2qi+CPwMrBUa31AKXUTsMFmUdVVDiZI3gyWfIh40d7RNChuTm64Obnh5+FH0p+SWHp4KYsOLOK3s7+x/MhycgpzuPfme0m8lMiKoyuI9I9kbMex+Hn42Tt0IRq8as+jUEo5AJ5a64u2CenG2WQeRbGf74GULTA6GRzkr9vasCdpD+5O7rg4urD44GKmr5tesgS6t6s33QK78cldnxDUJAhPZ09Zf0qI63RD8yiUUl8DkzH2g9gONFVKfaS1fq9mw6wH/PvD2eVw6Qh4ybLbtaFLiy4lz5/r/Rx/vOWP/Bz/Mz+d/InjacdZH7+etNw0sgqy+GrfV6w/uZ5ugd0YEzGGlk1bEuEXgcmh8d0pFaKmVPXWU4TW+qJS6mFgFfASsBNofInCp6fx9cKvkijsxEE5MDB8IAPDB5aUFZoLyS7MRmuNm5Mb/9n7H76M/RKAqBZRLB27FGeTMzvO7aBncE8CPAKk5SFEFVU1UTgppZyAe4F/aq0LlVINb+2PqvDtAcoRLmyF1uPtHY2wcjI54W3y5oNhHwBwKuMUe87v4UT6CXIKc8jMy6TIUsTDSx4mqyALJwcnfN19aevTljERY3ig4wM4m5xJzk6mlVcrXJ1c7VwjIeqOqiaKOUA8sAeIUUq1AupsH4VNmVzBpztckhnadVkr71a08m51RZnFYuHbB75l65mtJGUlkZSVxMELBzl04RAJFxPIzMvkjv/cgYNyIMAjgLY+bRkQNoChbYbSyb8TLo4uOJuc7VQjIeznuhcFVEo5Wld8rXNs2pkNkHUSCtLAK9IYLivqPbPFzPns8yw5tIQT6SeIz4hnd+Ju4jPjebXfq9x7873EZ8QzdfVUAjwDaObajN4hvXmhzws0dWlq7/CFqBHldWZXKVEopbyA14H+1qKfgbe01pk1GmUNsXmiMOdD5n5wDTTWgBINktaa5OxkTA4mXEwurD2xlvmx80nOTiYlJ4UT6SdQKFY8uILwZuGsjlvN2uNraebWDD8PPyJ8I+jSogvdArtJS0TUCze6euwXwH7gAevrR4EvgdE1E149Y3KBk/+F5I1w5257RyNsRClFgGdAyevRHUYzusPlH/mVR1eyIX4D7Zq3w0E5kJGXwbG0Y2TmZ5Kem45ZmwHYNnEbTiYnlh5aSnxmPMFNggnwCMDfw5+2zdvSI6iHdKyLOq2qLYpYrXVUZWV1hc1bFAD73oJ9r8Pwo9C0rW2vJeqd/KJ89ifv52DKQYa2GUqRpYiXf3qZJYeWkFWQVXKcr7svax5Zg5PJib/E/IWEiwl4uXoR2jSUXiG96N2yN22atZHhvaJW3Oitp63AC1rrzdbXfYD3tda9azzSGlAriSJtF6zuDt0/hPZTbXst0aDkFOaQeCmRpOwkMnIz6BTQCbPFzP9t+D/2Je/jYt5Fzl46S6GlkJ7BPfnkrk8AeOqHpyiyFNG2eVtGth9Jj6AeBDYJxNNZNoUSNeNGE0UX4N9A8QJH6cBjWuu9NRplDamVRGExw9JA8O0Nty2z7bVEo5NXmMf2c9vJzM+kR1APzNrMtB+ncSrjFAdSDpS0Sh6IeIDpfadjsVgY878xdPDtQKR/JB39O9Letz1tm7WluXtzubUlquSG+ii01nuALkqpptbXF5VSfwTqZKKoFQ4m8LsVkmNAW0AWqxM1yNXJlX6t+l1RtmjMIsC4rbX2xFqOpx0nzDuMAI8ALuZfJMIvgn3J+/jx+I9ojD8Ap/aayqOdHyUlJ4XnfnwOPw8/Wni0IMAzgEDPQO646Q46+BoTR7XWNHWVEVziWtXa4e6q9Z2eAz6s0Wjqm9Bx4OQFBeng0tze0YhGwsXRheHthl9RFkwwKx9eCcCl/EvsPb+XE+knaOPThsAmgWQXZBPgGUBKdgqHUw6TmpuKWZt5m7cxKRO7Enfx5PdP4uHkga+7Ly08WxDYJJCnop+ic0BnkrKS+O3cb3i5eOHj5kO4dzhh3mHSd9JI3Mg8ijNa65Y1HE+NqJVbTwDmAsjcJ8NkRb1j0RaSspJwMRmTCI+nH2fxwcUkZSWReCmRxKxEkrKS+Osdf6WTfydWx63m1Q2vXnEOZ5Mz/x71b7oHdufIhSNsPrMZd0d3PF088XXzxdfdl1tCbqGZWzM71VJU140Ojy1L41zCozSTM5g84My30Ob3xmsh6gEH5UBQk8t/3ES1iCKqRVSZx5otZlp5t+KutneRmZ9Jak4q8RnxxKXF0bpZawrMBfx4/Edm/Tbrms+ufGglLb1asv7ketadWIeXqxfert74uPng4+rDg50exNnkTEZeBs4Ozvh7+ONoupFfS8IWKvwXUUpdouyEoAA3m0RU32TsgZ3PglsghI6xdzRC1DiTgwlfd6OFUJ6Phn3EjEEzyC7IJjM/kwvZF0jNTaVrYFdyCnM4lHKI2KRYsgqyuJh/EbM24+TgRJ+WfVBK8ebPb7Li6AoUqiSRhDQNYd7IeTiZnFhzfA2pOan4e/jT3L05Pq4++Hn40dqnNSZlks56G6swUWitm9RWIPVW8D3g2ATi/yuJQjRaSincndxxd3LHz8OPNj5trnj/0+GfljzXWpOZl0labhohXiEUWYqYHD2Z3iG9Sc5O5kLOBVJyUnBQDuQV5ZFVkMXnuz7n17O/XnHOMO8wFt+/GIDn1z5PfEb8Fa2Vdr7tmNx9MiYHE6uOraLIUoSXqxdNXZrSzLUZwU2Dadm0pSSZKpA23o1ydDX20U5YAkV5xmshRLmUUni7eePt5g0YfR13tb2Lu9reVe5nYsbHkJyTTFJWEhdyLpCem46DcqClV0uKLEVEB0bj7uROem46Zy+e5UDyAZKzkxnTYQwWbeGtn9/i9MXTV5zz1pBbmXnnTBwdHHlq5VMAeLl40cS5CU1cmtAjqAejOozCpEysOLoCdyd3mrk2I8AjgBaeLfD18G00W/NKoqgJQXfDqa/h6CyIeMHe0QjR4Dg7OhPSNISQpiFlvv/ekIq3xtkyYQvpeelk5GWQmZdJel46TV2bEtQkiEJLId4u3py9dJakrCSyC7LJLsgmNSeV6KBotNZMXD4Ri7Zccc6xHcfy2m2vgYaJyyfi4+ZDc/fm+Lr74u3qTa/gXnQP6k6RuYjdSbsJ8QqhffP29XLdL0kUNaHlaDjcDRKWQvtpIJ1xQtQpIV4hhHiVnWQAfnj4h2vKtNaYtZkicxE7nthBZn4mGbkZJOckk5ydTBufNriYXEjNTaXAXMChC4dIz00nMz8Tjeap6Kfwdfcl8VIiIxaOAIxBBCFNQ2jn044nuj9Bv9B+pOels/bEWmPEmLMnHs4eeDh50LZ5W3zdfbFYLJi1mSYuTezWgpHfaDXB0RUGroLs01CYDiY/e0ckhLhBSikclSOODo50Dexa7nGtac3O3+8seV1kKeJS/qWSfptw73CWj1vOqcxTxKXFcST1CMdSj3E64zTnmp0jNimWP67+4zXn/dugv3FH+B1sPbOVZ1Y/g4NywNvVu2To8av9X6VrYFfi0+P5Kf4nTMpE98DuDAgbgJPJqUa/F5IoaoqrP+QlQ24iuEqiEKKxcnRwvGLuiI+7DyPajyjzWK01HXw70L9Vf7ILsskqzCKrIIucghw6B3TGz8P4XfJqv1fJzM8kJSeFlOwUUnNSyczP5OzFs6yKW8VbMW8BMLn7ZHq37F3jieK6J9zVZbU24e5qxz6Dnc9Aj4+h9cTav74QolGxaAv5Rflcyr9EkS7CzdENb1fv6x7JZYsJd+Jq4b+DuE9h+xRwbwmBQ+wdkRCiAXNQDrg5ueHmZNtpbY1jbFdtcXSB25Ybt6F+HgmJP9o7IiGEuGGSKGqaewgM3gSuAfDzPca+FUIIUY/Vi1tPSql7gbuBpsA8rfUa+0ZUCY9QGPorHJkJylGWIRdC1Gs2/+2llPpCKZWslNp/VfkwpdQRpVScUuqlis6htf5Oa/0EMBkYa8t4a4xbC+j4CuhCSN0O5iJ7RySEENelNv7MnQ8MK12glDIBHwN3AhHAg0qpCKVUJ6XU91c9/Et99FXr5+oHJ0/IPgVr+xiztoUQoh6y+a0nrXWMUirsquKeQJzW+gSAUmohcI/W+h1g+FXHooyxXjOAVVrr+nXTv8UQaNIe9rwM3pEQONjeEQkhRLXY68Z5MHCm1OsEa1l5ngEGAWOUUpPLOkAp9aRSaodSakdKSkrNRXqjnDzh9nXGMuSbRkPSRntHJIQQ1VIveli11jO11t211pO11rPLOWau1jpaax3t51fHZka7B8LANeDcDDYMhgu/2TsiIYSoMnuNejoLlN5GNcRa1nA1bQuDt8LpBWByA3M+mFzsHZUQQlTKXi2K7UBbpVS4UsoZGAcst1MstccjGNpMBsxwYr4xGkoIIeq42hgeuwDYCrRXSiUopSZqrYuAKcCPwCHgG631AVvHUic4eRrLe8S+DOsGwNlrlzcWQoi6RBYFtJfUHRBzj7HibKc3IfIVe0ckhGjkylsUsF50ZjdIzaNh8C/gfxvs/TPsnm7viIQQokySKOzJsxUMWAkho8DBCQoy7R2REEJco16s9dSgmZyh7//g4iHIPgmOHY2kIYQQdYS0KOoCBxO4BsHhj2Dv67IulBCiTpFEUVe4NIOCNDj4DqyKhIvH7B2REEIAkijqDqWg/3fQ6wvIOQtrekHcZ/aOSgghJFHUKUpB68eNjY/cQ2H385B5CCxme0cmhGjEpDO7LmoWBYO3GFupmnMgcx84+4JHiL0jE0I0QtKiqKucPCB0NDS9GQ7MgBVtYNtEKMq1d2RCiEZGWhR1naMHdHoDtBlOfGEMne1Z5gK6QghhE9KiqA+8boZ+/4PQsRA3B7Y/bezDXZAByZvtHZ0QooGTFkV90v1DYwht4mrI2A87noKULTA2DwrSwZwHnmH2jlII0cBIoqhP3FrAwB/h4lEwZxlJAsBSCEsDjecPNbxFHoUQ9iWJor5RCrzaQ27i5bLCi+Db+8oyIYSoIZIo6iu3QBhxFBLXQXY8XNhqlBdkgrOXXUMTQjQs0pldnzVpC20nw9GPL5ftf8t+8QghGiRJFPWdUnDLvy6/PvwP+8UihGiQJFE0BCZHGHnq8utjMs9CCFFzJFE0FJ6hMCIOPNvAnlcgfZ8xGqoyWsPXCg5/aPMQhRD1kySKhqRJa+j6N2NOxarOsHF45Z/R1gUHd//JtrEJIeotSRQNTcvRxl7cAElr4GJcJavPyrwLIUTFZHhsQ+TX25itnbQOijIhYy84ekLTttceqy3G145/rt0YhRD1hiSKhsrkAsF3Q1E2/PZ7SPoJhu2CS0fBpzs4eV4+TmZzCyEqILeeGjpHD2gz2ei3+C4YfhoAmx+48nbUvv8HiWvsFqIQom6TFkVj4N8X+i2BE/PgzBIIvQ8yYiE/HUxusO818O4MgUPsHakQog6SRNFYBN9lPAAKs4xbUJvHQGGmUZax136xCSHqNLn11Bg5eYJPN2Ml2psmXi7PPmu/mIQQdZYkisbMtxf0mgsOLhDxMvwYDXkX7B2VEKKOkUTR2CkHGJcHgUMhLwmOfQJF+faOSghRh0iiEAb//uARBvteh/95wOpoyDwEhz+Acz9W/nmLuZKJfUKI+ko6s4VBKRi0Ec4shaw4SFxr7HGx6znj/QfNRuujPBvvBOdm0HdRrYQrhKg9kijEZR6t4OY/Gs8LLoI5+/J752PAtyc4upf92YIMUPLjJERDJP+zRdmcmwJNYeRJiP8vODWBPX82Wg1NI8C/H7gFXD5em0GZ7BauEMJ26kWiUEp5AD8Db2itv7d3PI2KZxhEvgrZp+HUIsiz7svt2AR6/xtC7jFuW6XvMh5CiAbHpp3ZSqkvlFLJSqn9V5UPU0odUUrFKaVeqsKppgPf2CZKUSUeoXBvAow4AbetAGdv2DTKmOmdceDKY7PPSMe2EA2IrVsU84F/Av8uLlBKmYCPgcFAArBdKbUcMAHvXPX5CUAX4CDgauNYRWUcHKBJuPHw7QsnPoem7cFU6p/mawdAG+tLhd4PAQONFocQot6yaaLQWscopcKuKu4JxGmtTwAopRYC92it3wGu2WlHKTUA8AAigFyl1Eqti9fGFnbj4g0dnr/8esQxOP4F5J2HE19A3Gzj4R4CXf4KLQaDWwu7hSuEuH72mEcRDJwp9TrBWlYmrfWftdZ/BL4GPisvSSilnlRK7VBK7UhJSanJeEVVNGkDUX+FW+YZy5bfuQdMHpCfBlt/B1sfg/S9xtyMTfdByi/2jlgIUUX1ojMbQGs9v5L35wJzAaKjo2WDBXtr1hnGZoG50NhpT2twamrsj5G4xujbMLlDyEhjb+/uHxod5CYPMJXxY1mQCRcPG8uOCCFqlT1aFGeBlqVeh1jLRENkcjI2UAoZboyg8u5oDLnt9AaYc+DUQjjzLVw8Ykz0ixkJB94Bc8GV59k0GtbcAmZZXkSI2maPFsV2oK1SKhwjQYwDHrJDHMJeXH2h0+vGA4y1pRRgzjOSx55XYO/r4BZkHBv2CKT+ZhxrzjF25asJJ7+C2BfgnjPgIHNAhCiPrYfHLgC2Au2VUglKqYla6yJgCvAjcAj4Rmt9oKLziAbO0cX45e/sBbevh1sXQOuJRutDmYwE0mKwcWzm4Zq77vbfQ27ilTPQhRDXsPWopwfLKV8JrLTltUU95eAAYeOMR2mJayFhKRRlwaXjEPsiFOVCq/sheLQxk7yyYbjmAihIuzz6KuReiP8KZBCdEBWqN53ZopELHAyjL0BhOljyIWk9FGZA4iqM6TbWRQ09WoFrkNE3crUfIiDr+OUFDj3CjHJLUa1VQ4j6SBKFqD9cmxsPgPvTjRng+98yRlNZisC5udHyOPBXCB5pzCZ3CwAXf2je00gSAEU5xi5/B942XsuEQCEqJIlC1F8eLaHXZ1eW5ZwBkxscnQWUuqV024rLz2OnQ/jvLr/OSzEWO6xoGXUhGjFJFKJhCRpmPCyFRosjNxHyU8C7y+Vjjn0ChRfBLRhyz8IPHaDTm9DpNfvFLUQdprRueHPToqOj9Y4dO+wdhqiLcpON5GFyNkZUrYw0kkqzbnDnTntHJ4RdKaV2aq2jry6XFoVoXNz8jUex21bChsEQ/ijkJsH5DcaMcRcfoy/j7PfQ4U/GZMGalLDMWM6kY1UWTxbCvqRFIUT2GeNWlCUP1g+GgvSrDlBw8zRo9RD4dKuZzu+1/SFlk7EuVnXlJkL2KfC95cbjEKIUaVEIUR4P64oyFjPc8TPknYP8dCNhpG2HjP1w+AM4/A9w9AQ09F0MnjfByf9Adrx1578OcNNj5W8XW1rKJuPrkVnQrAv49696vEuDjK/Xk2SqK223MVosdIztryXqLEkUQhRzMEGzTkCnUoV/ML5kn4ZT/4OLBwAFusiYqJexB5I3Q9FFYzvYXdOMhQ67vm/sIX7yP4DFSDCO7sbCh56tL59+57PG16r+0s9LvvF6Vkf8V8Zy8ZIoGjVJFEJUhUcoRPzp2vLblhtftTb6MxKWgaOHMQNcm+HgO1CYeeVnAodde578NKMFk5MAzj5GK+fQ+xDxojFPpFj2qZqrU1Uc+9RYX0trmW/SiEmiEKImKAUhI4xHafcmQOElKLpk/XoRHFzBwRkubDWG7rYYBNknjUSz740rP+/gAm2ehJTNxvazyZtqq0YGc47xVZuNFpJolORfXghbcvI0HgReWX7bMuOrucC4haUARy+jn6Mg1diH/PQiY+mSnDOw9RFjccTSVrSH5r3g/Hoj4Yw8abRmck4bK++6NL98rKXIGAaszcYtsOpOLrQUwsVD4B5qLN4oGhVJFELYk8n58vOmbYxHsR7/NL6aC+GODZBzzugbOfeDMUor56yxAVTx6rffldooMnQsdHzZeP/EF8YtrWKd/wLtngXnJlWPUxfBys7GUihDf61+PUW9JolCiLrO5HTlUNhWD1x+ri3GqKzj88A1ADL2Gi2RokvGSCyPVtCs+5WJYu+rxgOMhRFvX2dMPjyzxPi8ydWYS5K88fJnLNaNpIr3BalI4UVY1dVYLt635/XWWtQhkiiEqM+Ug7HtbPRHpQoXXn7aeoLxKJaTBEk/GrsIXjpiJARzHmCBtB1wbqXx2lJqJ8HgEXAp7vLrgkxjsuC5H8C9pdHH4ugBLr5GPClbIesErOkFLYYY62yVbjmJekcm3AkhrqS10adhzjX6JvKTjZbLuZXgHgJN2xkTBos7usHou+i/xHj+y+/g4sErzzk62dhLxDO89upRE/LTYG1f6LcYvCLsHY3NyYQ7IUTVKGXc7ire06N4aXevCKOvwlIAPT4x5li4h4BXJ3B0A7cQQEPrSXDqK6MlUjyK66fbIXO/8bzr++B/m/UXrwPkJRkr/prcjfM42ODXUn4qpO+BFrdX73PZp4xO/JRfGkWiKI+0KIQQtpObYiQNzzYQU2rocOj9EDHdaLmsuWopEuUIN42Hdk8bQ4q3PmoME3awbplrcoGg4RD2KJz73uhbaXkvKGdj0mSLwcZclLzzkHvOeL7zWUhaC3fuMXZKbPsUuPpVHn/WSVh+E3T/CNo/W5PfmTqpvBaFJAohRO3JuwCZB4w5IR6twJxvtEzMedZbXfnG1+a9IGAAFGTAnleM4yx5xldznrFUSl6KMaHxan0WQpM2EL8QDr9fdhy+fcA92JqAnIzO/JB7oHkPSNoA51aAdyc4/7N1F0Xg/ovgVI2RYvWQJAohRMOitbHib36K0YluKTJujXmEGhMas04ZSakw3Wh1pO4wJjyCseshGLfRtNnog+nwvDFhcvU1vycNHZ6HiD+Di3etVM8eJFEIIUR5tAa0kTC0xWjVFGYafRRHPoIz314+tvPbEHw3pMfCnj8bLRIHZ+tXJ7j5T8Yqwxn74PgX1vecjZFfTk3hpgnQ9GZAGf0xJrc6szyKdGYLIUR5lALU5RnrJmdjBrpHKPj3g6I8Y+Ji9mnw7WXMUXELMiYgWgqM0WGWAuPh6GHc0jLnGcOELYXGLTVLoTHyy6+/kSBOf2u9daYu9784uEDv+cYosrPL4ORXxrHK0frVyegvcQ8yzmXyMAYAmNyM+S82IolCCCEq4+gK7Z66sswj1FhipTxNWht9KaVpizVxFELgEKPMnFvqkW/cFnN0M5Z0cfW3jjSz3lYzZxlJJ/ecsWjkqVJzZlxbQNunocM0I1nVIEkUQghRW5TD5ZFbAQOMR3naP208ymIxG4tFekXAiX9DVhw0bQ8ewUbHfA2TRCGEEPWNgwkCBhqPyP+z/eVsfgUhhBD1miQKIYQQFZJEIYQQokKSKIQQQlRIEoUQQogKSaIQQghRIUkUQgghKiSJQgghRIUa5KKASqkU4NR1ftwXuFCD4dQHUufGQercONxInVtpra/ZqKNBJooboZTaUdbqiQ2Z1LlxkDo3Draos9x6EkIIUSFJFEIIISokieJac+0dgB1InRsHqXPjUON1lj4KIYQQFZIWhRBCiApJohBCCFEhSRSlKKWGKaWOKKXilFIv2TueG6GU+kIplayU2l+qzEcptVYpdcz6tZm1XCmlZlrrvVcp1a3UZx6zHn9MKfVYWdeqC5RSLZVSG5RSB5VSB5RSU63lDbnOrkqp35RSe6x1ftNaHq6U+tVat0VKKWdruYv1dZz1/bBS53rZWn5EKTXUTlWqMqWUSSm1Wyn1vfV1g66zUipeKbVPKRWrlNphLau9n22ttTyMfhoTcBy4CXAG9gAR9o7rBurTH+gG7C9V9i7wkvX5S8DfrM/vAlYBCrgF+NVa7gOcsH5tZn3ezN51K6e+gUA36/MmwFEgooHXWQGe1udOwK/WunwDjLOWzwb+YH3+FDDb+nwcsMj6PML68+4ChFv/H5jsXb9K6v4c8DXwvfV1g64zEA/4XlVWaz/b0qK4rCcQp7U+obUuABYC99g5puumtY4B0q4qvgf4l/X5v4B7S5X/Wxu2Ad5KqUBgKLBWa52mtU4H1gLDbB78ddBaJ2qtd1mfXwIOAcE07DprrXWW9aWT9aGB24HF1vKr61z8vVgM3KGUUtbyhVrrfK31SSAO4/9DnaSUCgHuBj63vlY08DqXo9Z+tiVRXBYMnCn1OsFa1pAEaK0Trc+TgADr8/LqXi+/J9bbC10x/sJu0HW23oKJBZIx/uMfBzK01kXWQ0rHX1I36/uZQHPqWZ2BD4EXAYv1dXMafp01sEYptVMp9aS1rNZ+th2vN2pRv2mttVKqwY2NVkp5At8Cf9RaXzT+eDQ0xDprrc1AlFLKG1gK3GzfiGxLKTUcSNZa71RKDbBzOLWpr9b6rFLKH1irlDpc+k1b/2xLi+Kys0DLUq9DrGUNyXlrExTr12RreXl1r1ffE6WUE0aS+EprvcRa3KDrXExrnQFsAHpj3Goo/iOwdPwldbO+7wWkUr/q3AcYqZSKx7g9fDvwEQ27zmitz1q/JmP8QdCTWvzZlkRx2XagrXX0hDNGx9dyO8dU05YDxSMdHgOWlSr/nXW0xC1AprVJ+yMwRCnVzDqiYoi1rM6x3neeBxzSWv+j1FsNuc5+1pYESik3YDBG38wGYIz1sKvrXPy9GAOs10Yv53JgnHWEUDjQFvitVipRTVrrl7XWIVrrMIz/o+u11g/TgOuslPJQSjUpfo7xM7mf2vzZtndvfl16YIwWOIpxn/fP9o7nBuuyAEgECjHuRU7EuDf7E3AMWAf4WI9VwMfWeu8DokudZwJGR18c8Li961VBffti3MfdC8RaH3c18Dp3BnZb67wfeM1afhPGL7044H+Ai7Xc1fo6zvr+TaXO9Wfr9+IIcKe961bF+g/g8qinBltna932WB8Hin831ebPtizhIYQQokJy60kIIUSFJFEIIYSokCQKIYQQFZJEIYQQokKSKIQQQlRIEoUQZVBK/WL9GqaUeqiGz/1KWdcSoq6S4bFCVMC6TMTzWuvh1fiMo7687lBZ72dprT1rIDwhaoW0KIQog1KqeFXWGUA/6z4A06yL8L2nlNpuXev/99bjByilNimllgMHrWXfWRdxO1C8kJtSagbgZj3fV6WvZZ1J+55Sar8y9h4YW+rcG5VSi5VSh5VSX1lnoqOUmqGMPTj2KqXer83vkWg8ZFFAISr2EqVaFNZf+Jla6x5KKRdgi1JqjfXYbkCkNpatBpigtU6zLq+xXSn1rdb6JaXUFK11VBnXGg1EAV0AX+tnYqzvdQU6AueALUAfpdQhYBRws9ZaFy/nIURNkxaFENUzBGMdnViMZcybY6wTBPBbqSQB8KxSag+wDWMxtrZUrC+wQGtt1lqfB34GepQ6d4LW2oKxPEkYxpLZecA8pdRoIOcG6yZEmSRRCFE9CnhGax1lfYRrrYtbFNklBxl9G4OA3lrrLhhrMrnewHXzSz03A8X9ID0xNuQZDqy+gfMLUS5JFEJU7BLG1qrFfgT+YF3SHKVUO+uKnlfzAtK11jlKqZsxtqQsVlj8+atsAsZa+0H8MLazLXdFU+veG15a65XANIxbVkLUOOmjEKJiewGz9RbSfIy9D8KAXdYO5RQub0FZ2mpgsrUf4QjG7adic4G9Sqld2lgiu9hSjP0k9mCshPui1jrJmmjK0gRYppRyxWjpPHddNRSiEjI8VgghRIXk1pMQQogKSaIQQghRIUkUQgghKiSJQgghRIUkUQghhKiQJAohhBAVkkQhhBCiQv8fNsWyGfAxQjUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots()\n",
    "plt.semilogy(s1[:,0], alpha = 0.2, color = 'orange', label = 'Actual (Non-stiff)')\n",
    "plt.semilogy(s1[:,1], '--', color = 'orange', label = 'Predicted ')\n",
    "plt.semilogy(s2[:,0], alpha = 0.2, color = 'green', label = 'Actual (Stiff)')\n",
    "plt.semilogy(s2[:,1], '--', color = 'green', label = 'Predicted')\n",
    "# plt.yscale('symlog')\n",
    "plt.xlabel('iterations')\n",
    "plt.ylabel('Loss')\n",
    "plt.legend()\n",
    "# plt.show()\n",
    "# plt.savefig('Compare_actual_and_predicted_loss.png', dpi = 500)\n",
    "\n",
    "# Non-stiff case Error\n",
    "error_1 = np.linalg.norm((s1[:,0]-s1[:,1]),2)/np.linalg.norm(s1[:,0])\n",
    "print('Non-stiff case prediction error:' + str(error_1))\n",
    "\n",
    "# Stiff case Error\n",
    "error_2 = np.linalg.norm((s2[:,0]-s2[:,1]),2)/np.linalg.norm(s2[:,0])\n",
    "print('Stiff case prediction error:' + str(error_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step size plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002073563428166721\n",
      "0.00017585155024532433\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbkUlEQVR4nO3df5RV9Xnv8fdzfgxEo/kBZBkYI9iIC+THaCbRS2svKFW066JtTdCExCArXHsVVqxXS5dNqvbaaEOjITGxpCC5NhUNtUoMIYlerJKgEXQSEaSiJToRI5BiYgzMnHOe+8fe58yZYQbOMGfmsL/781prFnP2nLPPd8+MH5959nd/t7k7IiISlkyjByAiIvWncBcRCZDCXUQkQAp3EZEAKdxFRAKUa/QAAEaOHOljx45t9DBERBJl8+bNe9x9VG9fOyrCfezYsWzatKnRwxARSRQz+3lfX1NbRkQkQAp3EZEAKdxFRAJ0VPTcRSR8nZ2dtLe3s3///kYPJXGGDx9Oc3Mz+Xy+5tco3EVkSLS3t3PccccxduxYzKzRw0kMd2fv3r20t7czbty4ml+ntoyIDIn9+/czYsQIBXs/mRkjRozo9188CncRGTIK9iNzJN83hbuISIAU7iKSGmbGtddeW3m8ZMkSbrzxxrrse/v27UyfPp2WlhYmTJjAggULAGhra2Pt2rWV561Zs4Zbb70VgN27d3PmmWdy+umn88QTT/Dtb3+bCRMmMGPGjAGPR+EuIqkxbNgwHnjgAfbs2VP3fS9atIhrrrmGtrY2tm3bxsKFC4GDw3327NksXrwYgEcffZTJkyfz7LPPcvbZZ7N8+XK+8Y1vsH79+gGPR+EuIqmRy+VYsGABt99++0Ff27lzJ+eccw5Tpkzh3HPP5ZVXXgHg05/+NIsWLWLatGmcfPLJrF69utd979q1i+bm5srjyZMn09HRwec//3nuu+8+WlpauO+++1i5ciVXX301bW1tXH/99Tz00EO0tLRw0003sWHDBubPn89111038GMd8B5ERPrppu88z9bXfl3XfU4cfTx/8z9OO+zzrrrqKqZMmcL111/fbfvChQu5/PLLufzyy1mxYgWLFi3iwQcfBKLg3rBhAy+88AKzZ8/mkksuOWi/11xzDeeccw7Tpk3jvPPOY968ebz73e/m5ptvZtOmTXz1q18FYOXKlQC0tLQc9LX169ezZMkSWltbB/CdiKhyF5FUOf744/nUpz7F0qVLu23fuHEjH//4xwH45Cc/yYYNGypfu/jii8lkMkycOJFf/vKXve533rx5bNu2jY9+9KM89thjnHXWWRw4cGDwDuQwVLmLyJCrpcIeTJ/97Gc544wzmDdvXk3PHzZsWOVzdwfghhtu4Lvf/S4Q9dUBRo8ezRVXXMEVV1zBpEmT2LJlS30H3g+q3EUkdd773vfysY99jOXLl1e2TZs2jVWrVgHwrW99i7PPPvuQ+7jllltoa2urBPu6devo7OwE4PXXX2fv3r2MGTOG4447jt/85jeDcyCHoHAXkVS69tpru82a+cpXvsLdd9/NlClTuOeee/jyl7/cr/394Ac/YNKkSUydOpXzzz+fL37xi5xwwgnMmDGDrVu3Vk6oDhUr/4nRSK2trT6gm3UUOuChq2D6Yhjxe/UbmIjUzbZt25gwYUKjh5FYvX3/zGyzu/d69jWMyv3NV+G5++GVjY0eiYjIUSGMcPdS9G+p2NhxiIgcJcII93Kol0NeRCTlwgh3L3b/V0Qk5cII93LlXlLlLiICoYS7KncRkW7CCPeSTqiKSG1uueUWTjvtNKZMmUJLSwtPPfUUd9xxB2+//XblORdeeCH79u0DYOnSpUyYMIFPfOITHDhwgJkzZw75nPUjEcbyA6rcRaQGGzdu5OGHH+aZZ55h2LBh7Nmzh46ODubMmcPcuXM55phjALot0fu1r32NRx55hObmZp588kmga7mBo1kglXux+78iIr3YtWsXI0eOrKwVM3LkSFavXs1rr73GjBkzKjfJGDt2LHv27OHKK6/k5Zdf5oILLuC2225j7ty5PP3007S0tPDSSy818lAOS5W7iAy97y2G15+r7z5PmAwX3HrIp5x33nncfPPNjB8/npkzZzJnzhwWLVrEl770JdavX8/IkSO7Pf+uu+5i3bp1la+deeaZLFmyhIcffri+Yx8EgVXumi0jIn175zvfyebNm1m2bBmjRo1izpw5lfXVQ6PKXUSG3mEq7MGUzWaZPn0606dPZ/LkyXzzm99s2FgGUyCVu2bLiMjhbd++nRdffLHyuK2tjZNOOqlhy/IOJlXuIpIab731FgsXLmTfvn3kcjk++MEPsmzZMu69915mzZrF6NGj63Jz6qNBGOGu2TIiUoMPfehD/PjHPz5o+8KFC1m4cGHl8c6dO3v9vNzOSYIw2jKuhcNERKqFEe6lQvyvKncREQgm3FW5iyTB0XDntyQ6ku9bGOFeDnWdUBU5ag0fPpy9e/cq4PvJ3dm7dy/Dhw/v1+t0QlVEhkRzczPt7e3s3r270UNJnOHDh9Pc3Nyv14QR7poKKXLUy+fzjBs3rtHDSI1BCXczuxj4Y+B4YLm7/2Aw3qdCyw+IiHRTc8/dzFaY2RtmtqXH9llmtt3MdpjZYgB3f9DdPwNcCcyp75B7ocpdRKSb/pxQXQnMqt5gZlngTuACYCJwmZlNrHrKX8dfH1xx5b7vt/sH/a1ERJKg5nB398eBX/XY/BFgh7u/7O4dwCrgIovcBnzP3Z+p33D7GlzUjnli++u8daAw6G8nInK0G+hUyDHAq1WP2+NtC4GZwCVmdmVvLzSzBWa2ycw2DfjseVy5GyV+16HWjIjIoJxQdfelwNLDPGcZsAygtbV1YBNf4157lhLFkubQiogMtHL/BXBi1ePmeNvQKnWFe2dRM2ZERAYa7k8Dp5jZODNrAi4F1gx8WP0UV+4ZVe4iIkD/pkLeC2wETjWzdjOb7+4F4Grg+8A24H53f35whnoI8fz2LCUKmusuIlJ7z93dL+tj+1pgbd1GdCSqeu4FVe4iIoEsHFbqassUigp3EZEwwr3Sc3dV7iIihBLu5dkyVqKg2TIiIoGEe9VsGVXuIiKhhHvVPHf13EVEAgv3jKZCiogAoYS7q3IXEakWRrhX2jKaLSMiAqGEu6stIyJSLYxwr15+QG0ZEZFAwl1TIUVEugkj3LtNhVRbRkQkjHDXwmEiIt2EEe7lee5afkBEBAgl3FW5i4h0E0a4d7tZh8JdRCSIcI9uCAWG6zZ7IiKEEu5F3SBbRKRaEOFeqpoKqcpdRCSQcPeqVSE7dYWqiEgg4V6Meu5R5a62jIhIGOGuyl1EpJtAwr26cle4i4gEEu5R5Z4zLfkrIgKBhHt5+QGAYqF4iCeKiKRDEOHu1eFeUriLiAQR7tWVeymeOSMikmZBhLu7wl1EpFoQ4d6tcldbRkQkkHD3rhkypYIqdxGRMMJdlbuISDdhhHtVz728QqSISJoFEe7mqtxFRKoFEe6USpTc4k8V7iIiQYS7eZFOcgCU1JYREQkl3Et0xOHuqtxFREIJ9yKdZIGutd1FRNIskHAvVdoy1VerioikVRjhTpFCuedeUuUuIhJGuHuJTstHD3RCVUQknHAvEIV7SW0ZEZEAwt2dLEUKFrVl0GwZEZEQwj1aNKwQt2U0FVJEJIRwj8O8ZNluj0VE0iz54R732IvlE6oKdxGRAMK91D3c1ZYREQkh3OPKvZSJwt0oUSp5I0ckItJwyQ/3cs89E82WyeAUFO4iknLJD/d4tkzJmgDIUqJQKh3qFSIiwUt+uJe6t2UylFS5i0jqJT/c4567x+GepUShqHAXkXRLfrjHlbtnq8JdbRkRSbm6h7uZnWxmy81sdb333aselXtGlbuISG3hbmYrzOwNM9vSY/ssM9tuZjvMbDGAu7/s7vMHY7C9qlTuXSdUi+q5i0jK1Vq5rwRmVW8wsyxwJ3ABMBG4zMwm1nV0tYhny1BVuXcW1ZYRkXSrKdzd/XHgVz02fwTYEVfqHcAq4KI6j+/wylek5lS5i4iUDaTnPgZ4tepxOzDGzEaY2V3A6Wb2V3292MwWmNkmM9u0e/fuIx9Fef328glVK9GpnruIpFyu3jt0973AlTU8bxmwDKC1tfWI07hYLES3xs4OA6K2jCp3EUm7gVTuvwBOrHrcHG8bUoVCdM9Uq5oK2ampkCKScgMJ96eBU8xsnJk1AZcCa+ozrNpVwl09dxGRilqnQt4LbARONbN2M5vv7gXgauD7wDbgfnd/fvCG2rtCZ0c0xngqpGbLiIjU2HN398v62L4WWFvXEfVToRCdULV81cJhOqEqIimX+OUHisWoLZOpnFB1tWVEJPUSH+6dhU4ALKeLmEREyhIf7sX4hGouH1XuOqEqIhJCuJfbMrmuee6dCncRSbnEh3t5KmQmXz0VUm0ZEUm35Id7sdyWicNdyw+IiCQ/3EvFaCpkNpfHMS0/ICJCAOFebsvksnnIZON57mrLiEi6JT7cS5W2TBYsSxbXDbJFJPWCCfd8LgeZrG6zJyJCAOFengqZzebBMvENshXuIpJuAYR7dEI1n6+u3NVzF5F0S3y4lwrltkw+7rmrchcRSX64l7oqd8tkyVmJgi5iEpGUS364l2fL5HJgGXKm2TIiIgGEe1S5N+WjtkzOXLNlRCT1Eh/uXopvs5eJTqhmTVeoiogkPtzLbRkyWbAMeXOt5y4iqZf4cPf4hCqWVeUuIhJLfLiXZ8tElXuWHK5VIUUk9RIf7l7sXrnnTOu5i4gkP9xL1T33qC2jOzGJSNoFEO5FShiYQSZDDqeotoyIpFwg4R4fRly56wpVEUm7sMI9o/XcRUQggHCnVKRkPSp3tWVEJOUSH+7uRZxs9CC+zZ4uYhKRtEt8uB9UuesG2SIiyQ/3XbyP/xw2IXqQyWgqpIgIAYT7qvxs7jrx76MHlcpdbRkRSbfEh3tHoURTttyWyegG2SIiBBDuncUS+WzXVMiMpkKKiIQQ7k5TrvsJVd0gW0TSLvnhXuheuWcpqnIXkdRLfLgfKJbI5yx6kM2T9aJ67iKSeokOd/forkvDKpV7niwFVe4iknqJDvdiyXGnqy1Trtw1FVJEUi7R4d4RnzjNl0+oZnJkvaAlf0Uk9RId7p2FKMSbqit3CnSqcheRlEt0uB9UuWebospdPXcRSbkgwr0pG8+WidsynUXHXQEvIumV6HDvLMSVe1VbJhPfU1XFu4ikWbLDvVy557qmQmYoAq413UUk1RId7gcOqtxz0WOK6ruLSKolOtwrlXvVRUwAOQq6SlVEUi3h4R5Phcx19dwhqtx1IZOIpFnCw71nW6YpeqwlCEQk5RId7h2VnnvXVEiAnFaGFJGUS3a495wtU27LWFFruotIqiU63A95QlWVu4ikWKLDvaOPqZA5tKa7iKRbosO986BVIbtmy+giJhFJs0SHe0ex56qQ0WyZnC5iEpGUS3a4F3r03CtXqBY0z11EUi3R4d7VlilPhay6iEk9dxFJsWSH+0GVezxbxjRbRkTSLdnhXixhBtlML5W7wl1EUizX6AEMxGdnjueqcz6IWRzu3aZCqucuIumV6HDPZIxhmWzXhqrZMqrcRSTNEt2WOUilLaMlf0Uk3cIKd02FFBEBQgv38toypqmQIpJuYYV71c06dIWqiKRZWOFeWRWySKfaMiKSYmGFe7ZryV9V7iKSZkGGe7QqpMJdRNIrrHCvmgpZVFtGRFIssHCPLmjKmSp3EUm3ul+hambHAl8DOoDH3P1b9X6PQ7w5nslrbRkRSb2aKnczW2Fmb5jZlh7bZ5nZdjPbYWaL481/Cqx2988As+s83sOPNRuHu9aWEZEUq7UtsxKYVb3BzLLAncAFwETgMjObCDQDr8ZPK9ZnmP2QzdNkqtxFJN1qCnd3fxz4VY/NHwF2uPvL7t4BrAIuAtqJAv6Q+zezBWa2ycw27d69u/8j70tG4S4iMpATqmPoqtAhCvUxwAPAn5nZ14Hv9PVid1/m7q3u3jpq1KgBDKOHbJ58RssPiEi61f2Eqrv/FphX7/3WLJNnGEUtHCYiqTaQyv0XwIlVj5vjbY2VzZFXW0ZEUm4g4f40cIqZjTOzJuBSYE19hjUA5Z67ZsuISIrVOhXyXmAjcKqZtZvZfHcvAFcD3we2Afe7+/ODN9QaZfOq3EUk9Wrqubv7ZX1sXwusreuIBqoyz13hLiLpFdbyAwCVK1TVlhGR9Aov3LN53UNVRFIvvHDP5KLb7KnnLiIpFl64Z/PktHCYiKRceOGeKbdl1HMXkfQKL9xVuYuIBBzuqtxFJMXCC/dMnqxukC0iKRdeuGdz5Lyg2+yJSKqFF+6ZqC2jyl1E0iy8cM/myXqBTl2hKiIpFmC4N6nnLiKpF164Z3JkXQuHiUi6hRfucVtGC4eJSJqFF+6ZPBlKlAqFRo9ERKRhwgv3bLxEfUnhLiLpFV64Z/LRv6XOxo5DRKSBwgv3bBMApspdRFIswHAvt2VUuYtIeoUX7nFbJlMq4K7pkCKSTuGFezYKd92NSUTSLLxwjyv3Jl2lKiIpFl64xz33HEU6taa7iKRUeOEeV+55rQwpIikWXrjHUyFzaE13EUmvAMO9qy2j9WWkV+7wo6Ww7TuNHonIoMnVe4dmNgv4MpAF/sndb633exxS3JYZn2ln96/3M+LYYTTlwvt/WKq89QbseBR2b4Ph74IzLodjRx75/jZ+FX74OcjkYO6/wsnT6zZUkd50Fktsf/03tL26jx1vvMW4kccypfldTHj/8QzPZwflPa2ec8HNLAv8B/BHQDvwNHCZu2891OtaW1t906ZN9RnEr3fxu69P5x2/e523fRglDKvPnqVB3sEBMuZ0epa8FSl4hgM09Xs/RvS7fowd4If+Ycayi/fxX8zxv+NVTqj3sEUqOoqlSpt4WC7DgULUVchnjb+/ZAp/cnrzEe3XzDa7e2uvX6tzuP834EZ3Pz9+/FcA7v6FXp67AFgQPzwV2F6HIYwE9tRhP0mRpuNN07GCjjd09Trek9x9VG9fqHdbZgzwatXjduDM3p7o7suAZfV8czPb1Nf/xUKUpuNN07GCjjd0Q3G8akaLiASo3uH+C+DEqsfN8TYRERlC9Q73p4FTzGycmTUBlwJr6vweh1LXNk8CpOl403SsoOMN3aAfb11PqAKY2YXAHURTIVe4+y11fQMRETmsuoe7iIg0nk6oiogEKHHhbmazzGy7me0ws8W9fH2Ymd0Xf/0pMxvbgGHWTQ3H+xdmttXMfmZmj5rZSY0YZ70c7nirnvdnZuZmlujpc7Ucr5l9LP4ZP29m/zLUY6ynGn6fP2Bm683s2fh3+sJGjLMezGyFmb1hZlv6+LqZ2dL4e/EzMzujrgNw98R8EPXxXwJOBpqAnwITezznfwF3xZ9fCtzX6HEP8vHOAI6JP//z0I83ft5xwOPAk0Bro8c9yD/fU4BngffEj9/X6HEP8vEuA/48/nwisLPR4x7A8f4hcAawpY+vXwh8DzDgLOCper5/0ir3jwA73P1ld+8AVgEX9XjORcA3489XA+eaWVJXIDjs8br7end/O374JNH006Sq5ecL8LfAbcD+oRzcIKjleD8D3Onu/wXg7m8M8RjrqZbjdeD4+PN3Aa8N4fjqyt0fB351iKdcBPxfjzwJvNvM3l+v909auPd2BeyYvp7j7gXgTWDEkIyu/mo53mrziSqBpDrs8cZ/up7o7t8dyoENklp+vuOB8Wb2IzN7Ml6YL6lqOd4bgblm1g6sBRYOzdAaor//ffdL3VeFlMYws7lAK/DfGz2WwWJmGeBLwKcbPJShlCNqzUwn+qvscTOb7O77GjmoQXQZsNLd/yFeq+oeM5vk7lq/u5+SVrnXcgVs5TlmliP6027vkIyu/mq64tfMZgI3ALPd/cAQjW0wHO54jwMmAY+Z2U6iPuWaBJ9UreXn2w6scfdOd/9PolVXTxmi8dVbLcc7H7gfwN03AsOJFtkK0aBe0Z+0cK/lCtg1wOXx55cA/8/jsxcJdNjjNbPTgX8kCvYk92PhMMfr7m+6+0h3H+vuY4nOMcx29zqtFz3kavl9fpCoasfMRhK1aV4ewjHWUy3H+wpwLoCZTSAK991DOsqhswb4VDxr5izgTXffVbe9N/qM8hGcgb6QqHp5Cbgh3nYz0X/kEP0yfBvYAfwEOLnRYx7k430E+CXQFn+safSYB/N4ezz3MRI8W6bGn68RtaK2As8BlzZ6zIN8vBOBHxHNpGkDzmv0mAdwrPcCu4BOor/A5gNXAldW/WzvjL8Xz9X7d1lXqIqIBChpbRkREamBwl1EJEAKdxGRACncRUQCpHAXEQmQwl1EJEAKdxGRACnc5ahiZu83s1VmtsnM/sPM1sfbm81sTh32f66Z/fPAR3rQfivji9/jnnq/h0h/KNzlaHMP8G/u3uru44FF8fZzidbGHqipRFc+1lv1+KYSrcEu0jAKdzlqmFmWaB2Vfy9vc/fnzOwPiC7Bv8TM2szs5Hh9kofiCv8nZnZqvI974ztx/cTMfm5mf9zjbaYSXdrOIfbxgJn9HzN73MxeiRdmw8wmxNt+ZmbXmdmOeHu38QGnAyf0fH0Nx3+amT0S/8XyOTP7ipl9+Ei/n5JyjV5/QR/6qP4A1hGtlfOPwO/32D4p/jwPPAr8Xvz4QuDu+POtwBfiz/8A+EmP/bcBow6zjxeB/x1//ifA3URL7z4DnB5v/zrwYB/jawOuq359L8e5Fhhd9Xh4PPbTgHcAPwceaPTPQx/J/dB67nK0uQD4fWA2sM7MPunuDwKnAi/Ez7mYKAT/Nb7JVg54wsyGEwX3TfHztgLvKe/YzPLAu9x9t5l9tI99HEO0TPTt8cvywD7gT4GfuvuzVfuuXoXzVOCF+D1GAP/Q4/XduHvPe4POBJ519+fjsTZV7QMzW+HuV/T2DRPpjcJdjiru7sAGYIOZvQeYYmYbiJZDLcRPm0q0ouDy6tfG67q/6O7l2++dQdyCiU0AttWwj83uXow3TQG2xP+2VT11ElG1Xl6K9013L5jZFKL/CZR6vP5wWoj79GY2GnjL3X8UPz4GeNPMZgCzgL+pOkaRXqnnLkcNMzs/rlgxs/cRtVV+CIyl+700dwHnx3dmwswmx/fJnQp8wMyGm9mxRBX87VWvq/TbD7GPyXQP8SnAz4hu+DI+fm4LMLdqX9Xjq36P6tcfTgddt1j7AtENpMvOIOrjn+ruf6lgl1oo3OVocgmwzcx+CjwMfM6ju/G8AIw0sy1mNg1YQfS7uy0+gfmXccU/FXgAeIroxhBfL1e/seqZMn3to2e4TyKqvO8BWs3sOaJ1uXe6e/mmGZXxAfPoHubl13djZmvjCr3sX4A/NLPtRP9z2Ghmd8Rf+3B8TL/t+1sn0p3Wc5dgmNm/Awvcffsg7Pud7v5W/Pl1RL37v673+/Tx3t8A/ifwt8A6d39iKN5Xkk3hLsEws3bgAz4IN1M2s88R3Rauk+hOQX/hyb5frQRO4S4iEiD13EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRA/x+gcV/wyDSbjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots()\n",
    "\n",
    "sns.distplot(s1[:,3],hist=False,\n",
    "                 kde_kws={\"shade\": False},\n",
    "                 norm_hist=False, label='Non-Stiff')\n",
    "\n",
    "sns.distplot(s2[:,3],hist=False,\n",
    "                 kde_kws={\"shade\": False},\n",
    "                 norm_hist=False, label='Stiff')\n",
    "\n",
    "plt.yscale('symlog')\n",
    "ax.set_ylim([0,5e2])\n",
    "# plt.title('Step length vs iterations')\n",
    "plt.xlabel(r'$ Step \\, length: \\alpha_k $')\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig('prove_stiffness/Step_size_stiff_and_non_stiff.png', dpi = 500)\n",
    "\n",
    "# average step size\n",
    "\n",
    "# Non-stiff problem\n",
    "print(np.mean(s1[:,3]))\n",
    "\n",
    "# Stiff problem\n",
    "print(np.mean(s2[:,3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = np.loadtxt('prove_stiffness/Adam/non_stiff_Adam.txt', comments = \"#\", delimiter = \" \", unpack = False)\n",
    "s2 = np.loadtxt('prove_stiffness/Adam/stiff_Adam.txt', comments = \"#\", delimiter = \" \", unpack = False)\n",
    "\n",
    "s1 = s1**2\n",
    "s2 = s2**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26.386261790602333\n",
      "2183651.397699857\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABWsklEQVR4nO29a6wlyX0f9vv3ed/3vXPvvGd3dsndJYfkWCsv6HywDcGxpaWsmIItxFz7gyURIuiEThxDgBnISAzkgyEZAQLFhJUVQFAyHFK0IgtMsDKtPGx+MCVxRZHc9+7s7OzMndd9P877dHflQ3efU91d1V39PH3O7R+wO+dWd1dVV1f961//JzHGUKJEiRIlzga0aXegRIkSJUrkh5LolyhRosQZQkn0S5QoUeIMoST6JUqUKHGGUBL9EiVKlDhDqE67A0HY3Nxk169fn3Y3SpQoUWKm8Kd/+qd7jLEt0bVCE/3r16/j1VdfnXY3SpQoUWKmQEQfyq6V4p0SJUqUOEMoiX6JEiVKnCGURL9EiRIlzhBKol+iRIkSZwgl0S9RokSJM4SS6JcoUaLEGUJJ9EuUKFHiDCE3ok9ETxDR7xPRV4noy3m1W6LEmYChA6YpvhYnfLrsGcaAzr78uj6I3pYMpmG9V1oYtNOtb0aRiOjbBHyHiF73lL9IRO8Q0S2OwH8KwO8yxn4RwPNJ2i0xwzANoH/iL2cMGHbEz4z64vJhV1wuIzwyQtU9EF+TEQhjJL5fRnRl2HvPIqBeDNpAe0dQvwEc3RO38/g14OB9f/moBzz8AdA78l9r7wIP/kz8Lg9/ABze8Zd394Hju0BnV1zfzpvi73JwG9h9x1+uD60+iL79zlvWe4lweAfoHfrLhx1g75b4nfbfs/7zYtSz+jDq+a/1T4D+sb/cGFnvK4KoHgA4fWw958Xuu8DBB+JnMkBSTv9rAF7kC4ioAuArAD4D4AaAl4joBoA/AvB5Ivp/Afy7hO3mi2HXmhT6cFLGWDwOqmgwDTEREb1b/9hPgI0R8PhN99gAFiF98Gf+SX5w2yJOpuEuP3kA7L3rXzCnj4Hdt/yEpL0D7L1jEUhvuztv+jeWQdsiZINTd3nvCDj6EDh9JHiv14CTh+5y0wAevw6c3PeXP/qhvx7GrHE4eQAfhm2LgHqx/56/fsCqo7sH9A7815z6fGX2uA0EG+3pg0kfRRARVVN3/ytq3xBsuv1jYCTYDJzv0RVsfqaAQPJ9E21KR3eB4SmgSxgFUbmzIYqI+8H71pz14vBD4GTbP1/7J8Du2/7NfNSzxlvU51EH6B+J+5sBEhF9xth3AHhn4KcB3GKM3WaMDQF8A8BnAfwCgP+RMfZXAPz1JO3mDmdC8gtn712LiATh9JGfWHkX0qgHPPyRmzgyFp1r9MLL7TrcjJfoPfqR9S7efj/8gZ8wH9y2CDCP7oG1yLt77nLnPb2Lwll0XkLj3OclJiObAzQ8m4pDQGTl3sXtcJLe93fe0UtgnO/hJQTO/b5yu99dz3Jg9ncUccYzDUq3uiIwUFH6INr0gMm60yXcvndNTQFZyPSvALjH/b1tl/07AP8NEf0GgDuyh4noC0T0KhG9urtb4IXi5VpEp4HThxY36qCzb+30HY5AtncA5hF57L1ncY38PSKua9x2x71p9E8sbpcnQA5HLOJmvBPU2eQKMEELC2UCkTJxnDdQ2uNTjncYcgu4xhh7HcDPKdz3MoCXAeCFF14owPavCIfTHZwA1U3xPQ436eVOvRh55JvOUb+1bh+Te8Dyxcn1vXcBrQZc/KT1t8PljroANpRfIR4y/kSxOUDJ4k/KUcYlUtPiZJ2ThkjOzGKcJp251TsAVi6J6ysC1x4FqW88QZj+2GTB6d8HcI37+6pdNrtgNscrUsL4kPEEOrhtnSC8CJJ/Fgl5EQTmOaWELWyptYqMMEruF8mzg+6PAqcvUU5gzolUJE93EMaEuOrry59xxJ+dPf+1UKQ0L5x1EGWeORtiR6A8jwvfvCnOCSQLov89AM8Q0VNEVAfwOQDfyqCd/OCIVtqPgu8LheDD58plZABHhum1DnGUel7Z+ngD9RDH4am4HmdBenUUjpjKK/ZynhdtjCI4RMKnHGXifkZevCluck4fE89DD2Ty6bjwbrh5wnmXKIpRZwOLcvJxRKLeE5SzuYoU0wVBUpPNrwP4LoDniGibiD7PGNMBfAnAtwG8BeCbjLE3knd1xlC0I+54MQhk+l44i0BmAcHDucdLHJ0F5FWcjtuQnEy8HKQpUag69fs2lYgii8h25QX4rlHmlhJTkXZ9MZD2eslr/XlPUGlvoBkgkUyfMfaSpPwVAK8kqbtYIIQuDIeIjToAzmXdoWA4hLOzC6xetX7rAcdyGdJwtJEuvgIQT6DYMnrGJP2L0LZKPyNxuBEYARU475emUxeA1MdIFaFixIRWeSmgDMOgAhXC4MgzA491ORG61BaQSn9jEk2ZY5UMMpm0bxML8CQF/LbsBViEUqTxHVU4T5Hjlgxpj5ejn5CZOALxiLJXrxBUR1DbYfCOnczU2qEhUZiujFAS/WnDmYxFJD4qiy2Mq5MRa9lzIiciqzPhfQHk4+gsNm9/ijjuY+TEJNA0yUBWOq0oY5egD14DimnqMxRREv2kCIpRIgPPfTlOQyfb6fVJ2J8YxM3rcBWr3aIvgpSITpx5ELfOtFFr5dPOmUDxDTNKoq+CIIIZZWE6x3VeBBREFOM6R4mei0P0s+SCi6boVsX4ZJZwMyvS+xfY0qRE+ji7RF8fyAN88Qi9J8rijbjQYxNdQTtpRlrkEVc0ILVnLzGGlpPvZJ4bUCG9vFN8/2o9vboywtkl+jtv+mPOANZm4Njedg/E0Q55ZLlgguqWRZ4Mey5S+woLtNqMV3dRFr9WEZcXgRPPzYejAO8ahEytyNJGRp7gKSK3MAwzg503U6pI8SPPgF1vJiiKUxpJiP68QGWcC0SQhDAGQC0mc6GCSj3cqsbQgcp8kMuzyenLzAWVwizkiYIvxhIKSPANvRE780LQJlBbyK8feaGiIpKZn7V49oi+PnBHvuTx+HVxeSCmNBlU9BElMkBaojOVelJoq9JIXocLCU9ovvdWeEdR0p0SsXH2iL5IfNM7ykfGHPUYHXR/4U0hQ1B4sdaccHZph2FIiqEkLEcQvKKXaYijongZx0VOeq75EFIlgT4EDvNLVRYJxjBFWeacELEw1Bb9oalLFAdZWZFljcSiX5VwGPm859nj9L1IqrgVpagD0vmAsrqBSXTJIiyIEiVUEcsMOY05PgvrpCT66aF7ADx6XUIgEw506oGiOAR69QYdBWdhgs8bonrkKnyjMHNhJRTESspBHEewIjA2efQhKOdBijgbRP9424qRkYWHaZaTISiv6lTjpZxh5EmAVPU2SfvkfT5P35O0+z7LKMU7KSKXoFpZcFSeScDLFZNEBpxnRNWByBSduRL3FNuKU1eUrGtJPU7TEO/w7zhPRL8U76QJezAz0cDnGC9eNMFVJ/28OyE5kIYumCfiEIDUQx8XN+1falAKd5HDiaTk9DNArNydIZBlhioaZOEGVDFXHBWHeXsv7/sUxfN5jClZ7wTVkeoYFX8+nS2i78uDmgIc5UvhFleGKErcnCiQnXTS+m5RFbYuEUWW4seiz8viE8l5w/wTfVkmm9SRw+IKMuHkEYczipLJqgDZfyKjmpZnqmRsk2weSS3AUj+tTCmIoPyh1LsRvQsF6ENKmF/nrJ/4CevfpZTNoNoL/nrbC8Bid0L32wsAmcCirUNotwAQ0BwAVcNdz2IPIBZct/P3oA40hu5rYMBSz/38Qg/QPHWK2ubb6jUAoyJuH3C3020Cps0vtPpAxZTXm3d5YwjUdH95bQQ0Rv7y+hCoC+6vGEBr4C/XTGChr17umgd8ObO+PQB0G4DJnUSivjP/bTotgHEbUFUHmkPBM4r1OXPXgWweBdYXcE32XeL2T/YMP96dJsA4flflnQY1YFQLb0dWl2x+8HX164Buk+T/8B+QFeaf0z+LMGN8VorJybCiiw9mAGkOYcUjekv6fTTPSdmYA5KhTYNrV2iznk9okvnl9J2d8sGfpVvv5ef99V5+3nL+ckzfLj9vhXfYecP6+9KPWcf/g9tA/9hdz+M3JwlFZHU7f68+ARzfdV8zDeDRj6y/tz5mpb4TtXN8H+jsuMv4tpprwMZT4vYBS0z26IfW7+XLwPIF6/fhhxNdiajevMtl79nenaSk5MtPHgDtx/7y/rE1jt7yUQ/YfVu9XB9MvL75ckMHHr9m/T73UaCxHP+d+W+zeg1Y3Jxc6x1NwozwzwTWx82p9etAa31yrXsAHH0Yrb6ga+0d4OR+/P5d/HOApoU/w4/32pPAwsbk2uGdiYe7rA9UAS7dDG/n4AOgf+QvH3Ym+TtU28gIuRF9ItIA/E8AVgC8yhj7rbzaLgyylgsGhsQV5EHl748StbPQycQzgGm6CYsIKjlyR33Oj2B+ZMSJkNg5y4SawCJpO3kEZMzHQCLRWY2IvkpEO0T0uqf8RSJ6h4huEdGX7eLPArgKYAQg4yzgRUAWdvohxNZRJsZZSM3V6M/MMzJJdH7GNss8MKthHaaIpAK6rwF4kS8gogqArwD4DIAbAF4iohsAngPwnxhj/wjA30/YbvEQZL2RVnIWR2zBw8VJ2oqqOIkuojiunSXz1DjInaik7KFadI9XV58K2L8w6NO1fktE9Blj3wHgNX7/NIBbjLHbjLEhgG/A4vK3AdiCM0jPMUT0BSJ6lYhe3d0NiD0zS0grZEKYfbxDuOOYJ6qagwLFJARj5OghLROJycqTOvKpeI6mLSLIOptcXsnfpw3XRjpdP5csVPFXANzj/t62y34PwE8R0f8K4DuyhxljLzPGXmCMvbC1tZVB96aAQG/ACJ9gqkHWMrBPLxoiJ7mRiWtkoqEcFvvpQ2+j3M8YG183Ay92HmG6kjAUjgEpTrgFGXLbZhljXQCfz6u9/OEhfqoftrkKtBVPAmGBsaLKjONOPhcnO0dE34UEIgSZ+KFwBCptEPIXtyi2V6ixn25fsmAd7wO4xv191S4rkTWCQjGnOen5lHdpJcquL4ffkxZUlLS82MGl7/BY44jKVURMhSJCAYjifd3I8RsmQk5jn1s0gGjIguh/D8AzRPQUEdUBfA7AtzJop1iYmphDNcrmvHLkHiQipvyzCuMlS3phSJxskloEZWFRJK4w5foUMTUFZ1ZrQzKOU45dldRk8+sAvgvgOSLaJqLPM8Z0AF8C8G0AbwH4JmPsjeRdjYFuBgHW4mCaWa6mxW1U0op1kxIiE0wVrlxWZ8Qx98nhM0LW0SoT1606bgn74FW015eS1ecCt4FIrYxmWKbPGHtJUv4KgFeS1J0KHM/BPOBTsnIf9vgusHguv76EmbQZOlCJ+OmjLvbGMtDNMJVkEpg6UKmF3+egdwg0BIQhikMbAOlizzKAHf/d9B5QX8yurTiYRuas3gGw/mSydmVork48co1h9HWWA+YgkEZRINvhJcjtKCs4upqO+CENjoOrgzfvK3I6x6hmiDKrGxVinRVnHOcEZ6Yc2yV1MUWMsRL5ruTRrgoKao5a4JVZUMgWW3Ml/FmXXD2liRZKVHI8Sk5LVqlCWF1JZCT3y0Ic8+/Ft8Urn1XssBM7FXHPpE3A4yDtRN79k3Trc2FGFOc5oCT6USH1XI3I6YsgUwAmgbMwc0vaERcJ7JtdBFChHpmTlIyIDSTEaCipRyajlyXxibVZlkSsRDyURD8yUrILdhZ6Ji7lIfVElkUnRGSnp5QUpNL6FerJ1aQyRltx5k1SK6EobWVpQJDlt4mSTEgElb7NsvXOmYSSvF4hjs0opdAMuSFJoLG8RFkJ75+lNJBxNsA8N/swR8I8EGfexem3dIOTWICdTNdtaX6JfmYKNJVJodK24J5Yys8Y7+lsSpma8EWsW+Xdo1q5ODHSvZBtykHObbFRIDFM0tg/WSOtGFUOcjutSdqRRQAtOf15QYEWtxBc/5wEK2kjSb3VZvg9/IZbqU9+q3CwPAGQ6k5ibtbKz0r6k+UzWSJtPxjZJh0bih7QWZmNpsUQpYz5JfrTXCD6AKEfNi/b7DwR9WgctZ/8xsBnc5Itbln9slPFLIl34iBtTjoooXsRNqhpc/qFNJiYZ6KfFVRk8Sr39I78ZVlw4CJuLMpkFCqcU4KLyDJJOYeoMumBZDxlVjcF8pqMhDzDMKgqkDMVJWX4beKsQdmYZB2hNCbmmOhnNDGSHEGFi1OB2CWByJ47SojfNOzBDQlHKOO4RxLiHrUv88y5R0l6M40+qH6rIpwIeGS2WRXnPeeY6IeglqU7ehLZsPe2kPtSCzCWch2ijF6Jm1KU0ToogkdkVkRNlfkoRKycGcI8vpMH80v0wyb71rPA0sX06vWWF4WDcTiXafYnLVNBmR7E5ckpEz9kbDYaeXwTeuSmTZyKEFJhaojITMw45pfoq2DpQoyHUpoUWSpy0yR2Ufopk4dK5aQR+ybzjE07HMAsIG3iNHN+IwoYyJiBs405JvoBH7lhx8nRNGD9uvX38iW1alVklVFM2YoeHmFsoaGwaGQbhIpVjwoRc1nsZIzCc3wF798gQs5lGTRJJFTVbxOHsYqT0pLPL134eTPXRF8RrXXg3EeA5YvAwmaCinjiHfOo3H6crN1Zg0wRKBU1REx2UdDMRZGgFBE05TmQRhgGmTI+WkeiP9La4P7IKXGQykkztwQ44Zhfou8MpigNn4yjXkxC9H2NxHhE0K9MFXE5eeRWW+JbRhKi304jXC7cpplSm/2sCUNG309ZpFV0piCgf1SRXysRG/NL9B1EWdQ1CXEKQ2oWNFkqG6fprBZVXpxjnB/es1etgYjlJVJBnAi0quaX00jkMkXMMdG3P4SI6Ket+Bt1s/vwmSp8I6AIE1umGyiCzXpUxBnP5lrq3Rgjz/j8sUJQxJG1K5z00kBBnbBkmGOibyOKeCcIWx+f/M466xVvXx5XP5AWMglCJoLCopQFsErU7Axx7tUM8w6n4Q2uSliVzUPjfJtpf7dptx+O+SX6LIDTD0JtQVLOxX3Ze8d9TbRgksQujxJtM2sOPCi+SonskURJGBsZBiebZWThCFkqcjOAQ0BXrqrdv3JFfu3y85ZDl5dTGXaQ2Q4fZPoWRc6ZmZw/wWSWiRSUFkIezlApIOmiziuX8rREZMrjk5MlThDmZMPLlegT0SIRvUpEP5N9axynf/l5YGlL7bHGUvD1lUuWiSePYVttQsRK6hBA2GUxbVIlcBlOdJmJqvS9VJDxwiyQ6Z0SVPuVp4Nbou9bIikSEX0i+ioR7RDR657yF4noHSK6RURf5i79YwDfTNJmjF5Ofq49af174VPBj1x+Pvh6Yxm49GPuMhXZe+bmgVkiA6I2bX3FVFHQTSIJago5EYB4SuNU7P5lmMNvEYCknP7XALzIFxBRBcBXAHwGwA0ALxHRDSL6awDeBJCSEXYIRBzOwoZF0CsKgbjOfTT4unOCcOALgCVoX2gSGjbhAq6rxO2fFUyNU54hzr2IfYqDOPqu9q64XLnunMZuBr5RojCEjLHvENF1T/GnAdxijN0GACL6BoDPAlgCsAhrI+gR0SuM+c1oiOgLAL4AAE888USS7jkVxnuusWx564ZFM7z8vGUPvH8rXjsi8BMnyEFlcJqtRYcKpr3ISiQAYSq6kd4h0FrLv91pQnUz0IdANarvSDRkEXv2CoB73N/bAP4CY+xLAEBEPw9gT0TwAYAx9jKAlwHghRdeSOKlE/9RB+vXrf/C4Ih79AGw+5ZVdnA7eftAvE1LVWkrm4iVxuzKXQsnc0+qyFVRsMZso7EsDmKX+VjFcEgsQlC9UTdc5xcZ+c/L3AOOM8a+lm+LOcnRiYJlmgNFZa8XgT4FGU2YqoDoF+3YWjginlF/unvA2rUEFeQ4TqrfRGYWXXRETQdaUGRhvXMfAD9Lr9pl+WJaROHy88DFm/4IgfvviSdNWD+DgrD1DotHjNNGc3XaPcgGRfhu0+pDHEcwWcTNNJDEp2YGkQXR/x6AZ4joKSKqA/gcgG9l0I4apmExo1WAi5/0WwEJ9QMFipNTxEkdNehWFMc2oJjvLIWqx2vBo4sqW+9w7xvEZWdlLh0LxZ9PSU02vw7guwCeI6JtIvo8Y0wH8CUA3wbwFoBvMsbeSN7VqCjI4IeZf0YlUoVDWjL0lL5X5ABqc4i08zKkTTAT98/TH9Vk7SUAJLfeeUlS/gqAV5LUnR4KYBt/+XnLi7f9GOjsTezTTx9l1KCq0vYsLZA5Pj1lHYaBmZMTVxptpR3cjZmY2eAC3vHUe5lb78zoSCmgaMd2rQKsXAYu3ZyUnT5UT3AtxZSOrXFC3aaKKZmKZmUdlGaOWl4Ukus6iGOokEL/Ztrp0YNh9lZK80v0HRRxQlx+HqgviZVTpplvmFseTnC1JHFzwiCKehqIgm3eqcHzXmkSfeVw3PM6thnh6MNp9yAV5G6ymR8KPqE3n7H+Pd62js5tW9SzMwX1hwO9p+5K76BoJ6rIKJpdfwDOmJXJ1JWvsUysiz/288/pF0GmH4TVq1YQNwdpcPlJnbPUGkl43Xt78ReLEFHFPSpisaRpAtNW5BYN3rHl/1b1EE8qRonqa1Cg+T2/RL9Ag6yEy88Dy5fC7/PC+55JzPWixDRRTUXnryjmcwmR1XwIrVegFBWBF0NqCYl+Uj2RT7nIewXP2LqSIWmwPy2hkGQ8xvmP5/wSfQdFlOnLsHzRIv6tDeD8jRjyb2QcjZDDtLnJwoRbmJJpalCdaTu0yRLY+7oR591SGI84uQCmzRTK2o/NTKljbon+9mEX7+0EJCApMtaftEIhbH40uidiKpM5Sh1TspbJHIqK1qibnyFxMhqmOFdditwIp7dZxTAOo1PQd8+BaZtbon/ccxZXupz+yDDRGwYfDfsjAwM9JWuMzWetaJ9FQ18QqGsWoawclRB3GZc56gXf7xOhcLGOkkYu7R2p3Vd4KL5vrKCE3PdM9dRY0M2Ew9wSfbLHniX4CKbJ8Nr2MQ47E87p3cenuOU5QXSHOkbGZBK997iNdx+lxLlV6yGRPj3vN+r5ywA15W6USR4W9XBaHrklkHwsi/4tUrCqkZ24MoPACWtKmFui7+DxiXqO0c5Ax2vbx+iPLC5dN60P9fh0ws2J9KTv73R8G0HqOH8DjAGnCyG5fpOICU4fBlwMWGhFSuQub1jxtpStjnyWJjnoQmZZbJNx30cGs9b1LI9RQswx0bc+Ks+Bi2CYk4/viITag+hmk7qR8SSqNvBo6WO4066HipdUwbiJP9BNGLoa9zPQTby3c2qNU+GIeEb9cTZTH8c4FJfLoEvyFHT34/VL1HZQcLI4TmAqohDrYoy6850/d/Y7uLOXk7GDCqawfuaY6FsgTS7vO+6O8OaDE3SH0w4pIEd/ZGD31CIUQ91afHqQXbLiJHp03McbD07Ht9896GL7qOurQzf99TmbTup6E+GJoSAmnlK5b1TOXWa/r34ijdZc0Hsoju3Jdrp9CEKcTSmiiWoheHzZ5p8D5pjoW5+Wn28D3cDD44ksrW0T+24I55zFZizi1u/ud7HfHoAxhvcen+K4N8KtnTYeHYebpPVGBgb2piCV33NOP/sda9LxNN3ZVPg+frDXwXHXIu66yTKN2muYLHCsDzpDvP0obQWymqy3O9R945NquwHXGLN0S0dd98Zw0reYFpNZ3+XeYQ/90Yw7ZhkqCm0PeP1SqsrXFE4usk1sipnA5pjoO5hwoXf3u9g7HY5l9l7IPnHaopuDzhC3dto46VsExiEmx70RHhz1YZgM/ZGJ7cPueA7xohjU/Snbtg97uHvQBf8WvaGBD/Y7MAz/+5IxBIVM6r7dr87QEuN8sNfB/eM0FFCCkwFjuL3XxaMT+Qa33xlipEf8FqZkI5TtXhLRyJ29Lt555LahHugmhn1r8b758AR7bY4o94/E9UdMIMLY5LTlnPgcPD625srQMDEwDNep0IvxiU1g6dIdGhjxu3/R5d1p9C9tBzYZZB72Uxzj+SX6AYPqXPJOf9Oe+IZApJEG+iMDh53JpjMYmTjoDPHOo9NIIiZWWxz/Phno8KotRqaJ/sjEXmcI3WDoDUSEIN47yjbMSbURPVQ94PUpI8OcnF4CcNzTxxuoanuPj07w4f5EttsfmXhvp+36DrrJsH3Yw1CiF7p70MWH9x+M/z7kOXGZuCYgzMZxbzSegwBw0B3i1m4bpmdMB7pgXAKGtWuf2NpDHZ2B7tNZ3T/q4f5hzM1c597T08+hbo5PiV44olWfdZ1g/vRGBg4l9ajguDvCsXd+eKxndk8HgQxHEE4C5t+9gy46tu5rZLJE75EW5pbotxqWKKNSmbximDnvkf1Bdk6ykbe997iNbc/i6tgLcBDhWM6aGwAsS4R7+x085ifrsI0P9jq4dzg5PhJgLaZMvJM99driEd0w8fC4L9QJqGBkmHj38Snu7HXw6KTvd7TjiMPOaR8f7rn1EacDHe/vtl1ElMfuyQAnvQnxc4h9uz/Ca9vH+GCvjdO+jt7IcN0HMDw87uO0PyljzG0N4mwSpslw/7DrGwMmIGz9kYm7+1084E5STrteJuTuQdc+1dn18dcFn9jZqPtDE7d3O/hgtwPGGB6f9jEyxAYPnaGOg45o4/IS6QkT8PajU+y2BzBNa27e3mu7+snjwXEPhsn880PgaLV92MNe21qT9496vjUUhrsHXWxL+mGYDKbJcNQbjb/paNCTzpv+yPBde8zPPw6maW3kH9iK40fHfey1B9aG3Tuw7/K3Y5gsUz3j3BJ958XcXJK1IpLY7mcFET0mbgUzNiECe/ae5LzbyCOqMFXEUSEJ1+PI7ntDAyPbAui4N0J7oI83UgA47Axxe6+Nvm5aIiMJjrojvP3w1OoDg4vA8sTVNDkpjed9DjpDmAwYDqNxbzSwxC/9oeGeJ1z97YEu5QqPuiN8uN9Fr9vGYXeI456OfZtg3dm3xFev3z/B7d3JJmaYgG7XP4ohSgzeyyf1EQHayBJRdYcGjjtD6Xs8OOpjX0j0g3HUHeHRSR939jsY2YwM48bu3kHPLYLyGg1wyvy+bvo444P2kDMg8IOBuedLCN58cILb++6N5p2Hx7iz7998TNNi3HiGSgTdZOMNU+M8bPnNngUssLsHXby/08lM4jC3oZWdAWbMvyLSEKeNDBO1itqeedof4Y6AE4gLXsdQlcmOeXi8Q0cmgzbqwNBa4tsNE3f22qhWop0Mto96GGltfOzqwpjU8ARp+7CHhs5w78ByIHvmvKWb6A51NBjgbMq9kQHIIjxzBOT9vbav3GAM2wfdidI1yHwRABizzFVDJkW1twujdS64LgB92xNbH/RAHmo8MkyLo16yFfn28N/ea1tvvmn93RmM8OBQYFbo6SMZAwD+MB26yfDW9jGubbSw5rlGacbtl8DZ0MkcglUarmt93bDGaMX9THuo4+FRH5cXdDjCy3sH0U0rDztD7J/0ATSx3FQjb7rO4A1x1xkYvvF2/moPdECa3Iph52SAzlBHo6YJmSvGgFu7bdQaHVxf8/fRMfKwaFj6p/O5JfqOkpL/bjKOKM4m4JWzBkF8TE4obVm9Bux/CG10AlQlm4/TRXNkuygTTvsWl2ourk7clj1wOAxnc+kMdKDJcSme+ztDA4en1sbCRgO88eAE9fbAdYw86YnHQDdMvL/TwcYg6abo9NVwW9kw94/9ztD6HmuTW2QiCOWWU2Ai+Co6A0XCzBGUyWsya2OpAHvtIdaaYnF/cJfFV7tDHe3TAc7bH/akp2P74QkaLQODkSkVGzAWTrocQtcbjlDVDTSqQZFG5b3XbQ46mljRfW9leAK9KmaIAOd7W8zCkLOY0w0Ttx6eomJvemE9sE4kFgl2DDiubUzazUoeMbfinTG9C7jmmGqqmEQWCQwMWNwMvEOGXpgiVgCZGaCjK3hw1B/XqwncyxkDPtx1LF/cfXPWZm9kQpWr0SUKXsNkGHjfj9tZ5XLqCSqSkxNJ7fTF4/nopG8rR5MvXWnTgrK+Z1weHPVcYUQqo+Aojrph4o0HE7NY08R4rN/f6eCAs1DabfdhmkB3YMAwWag1mPAdPI90Dx/h3Uftsa5LBYddSw+jBzlihu3MfA4DxV387oHb2kym8LcrDexXZ6iPT4lZY26JvjOYfon+BFnJzMKQtNWRzvD6/WM8HAQnUHYmkcjipjIUL37TNgHksbYgjvQpt5hxQ7aGHEWXdZP7GgU4K/V1S/zEQxt1cXu3PfY/GINTDD44mizQylBi789b1/B9Yk4/ZV9P/lW9Yp4gWBuaePF7q6GxT4G87eOePr7KGKC5RH3+5zoe/5HHp33cPehyyku3L4Vokx/XziVTZwxCyxXvOw2Gjm5DXal0as/DztAIty6Tgu+Ie1zaQx1DwxzreLKO1p51/bkRfSL6WSL6TSL6HSL6yazbG3P6zF+WJrEXWWKoPiP7tpO+B5udnjYuAwBGHHdnKU7dz43FBdxski3WuwddPDpSs44Y28xz9ZJHhqkbpit2EY+9zmByypLGqPGPgYijI8aEJ5KhYVnFvP3QvcmRIe7ThOhN1LhhZqNpsg79kel7D1n9miwML/P+YS3zw250xaxjRSJSbBII2kgt1tNxbzS2wHGDuf5Jgrv7XfRTCFHC68kGuomHR318uK8qAuTEoDHe6fFJP/OQLomIPhF9lYh2iOh1T/mLRPQOEd0ioi8DAGPs9xljvwTgiwD+dpJ2lcAcy4FJUdcmfjIHlsy64rVKU/Ts5J8TKhq1Koz6qmu97IoWFmdS6N2kvKZhsrhDTDHcQKW/N2kTwEnfsgsXce66wfULDMzePKx/AjjnCMqt+wddibWHbcnleS+eU3S+Ez9GcS2/HvpEiNaJSpjzQfk0wetZ/M8IYzTZJxYRs3LSG+LuflfYPq9nYMw6pRmMBXL6PEzJuMleVdUC57SvB27KUUOcM4bxGPV1w6Xv4U8f3vGTMRHV/mFom3xEgJMIlkdxkZTT/xqAF/kCIqoA+AqAzwC4AeAlIrrB3fJP7Ou5QKjE4ghN7Hoz3IxFdR9JnDpGy09gsP5cYAWaIbKXJwBMSY54fNrGWxyn7B43tYHQDMlGy6J5LJIxEB+RJPVourNoPWMSJ9uSBGGnvf7IEG6mUXwzAGuT/JC3z7eb9Z4I405Nka25CAazbPABD7GTiOj6uoH9dsApQxBz6Uigexl6woy0h3LTWbKvv3ffCWQXMCoSK5mRvTYcCyTeP+DWrl88Ko5+Ehbw0cR9xZN1WkhE9Blj3wFw4Cn+NIBbjLHbjLEhgG8A+CxZ+FUAf8AY+76sTiL6AhG9SkSv7u7uJumb9W/APUmOUSrybJmDh1McplTkId2giHxmcV6M2ru4tWPJuyecRLR3l9G1tx+dCs1iBR3l2rV+B4nZZApHTe9OlLWc8q0yEIc3qHDiB6ZNdCBk2NZGEaeASIQUGMHANgkVQdXChN9U+FOi2d4BAHyw2/U4kKUHURc/3O8ozR5n7P2nnAmGuokP9zkHtoDctd4wGA+P5PWe9kfWdYmJ6kA35eFYtOD1BITScgBA1+VdboxP67zM/q2HcsW6ijd6HGQh078C4B7397Zd9g8A/FUAP0dEX5Q9zBh7mTH2AmPsha2trYRdoUimlb6+KEztoOqdj+wN6OYQu/7IDNH4pwNrc2PYb/sJQ39oxHIIigwBJ86fMqwNklsNAQM7NEyAuYn4eCVF5HqDzDVFXpF3dv3imHcencqtewIglnHHAYscDpyN/8cXqtURQ4slvbLfGWCom1z/E64H+zt4rZi8uHvQxXuP/d+SMYBM67tMFP9upkboVyX4/l1uU7n1uB2ZyYxiwRQFuSlyGWO/zhj784yxLzLGfiOHBvl/3JdSqD6K118QR+voGXi5cxyx03DlqfFviwCFfFqmwxLvmLjnI3wMZm0S1I28ctsYA8hby5gBNtB++BsTflMt2JJJBnlUSibkcsMQaDYYE1G7oTJ/HHELbwRAAXGBJpWr98aZN5rAHMUcr0+WqbkKhTnnebDfnWzEMoZM6u8D+MaHzJE/1t+UA9plQfTvA7jG/X3VLpsCxJx+VFmqCF1VBxoBRIrcsON52HWzvjz+rY1OQkU+PERKYsaJTix5vGSiej1v7OcIDKwiIO4MYBVJsneuHplJZXVwJFx0lYFXyqiKdBdgOkdyz2nFMZKK0VfyEVXrd9tmWsSeyOmMSWVoidxEG6HDLfeH8cbrQNESyTFNlvghjuEYEYh1Z2qb0ge7bgsni7P3v9/eqVrfs9obsiD63wPwDBE9RUR1AJ8D8K0M2gkEC+D080KaO7pKtqz+5k0Y9VUAlgu8D965K+jfmD7E5r64eEGausO3l6CRRPFLej/mNxUr67KHRCEd4R1GozQjM2Y5BuIxdk5T/OkxKQIVwxK8LxDNWZB9owCLIy9HrygaAyZxloLaBrI7ACU12fw6gO8CeI6Itono84wxHcCXAHwbwFsAvskYeyN5V2N1EIzFt8tXIS5Bt8TlYpLsFaPlazDqq9BbHo9dxUqdIGexOyS1ohETcZmTGJOKgNiEG8vaiyUGsuAxHoXlMIikx1AjOHEh/25AtA1HJMOT+26o1BPbPUdhnlV7YSdNJvwZhFY9KBRFfCSKvcMYe0lS/gqAV5LUnRS8XLM3MrDUmNswQ26QhtHKkwCAavdxJk3EY7R5bpVz5pLIXEmXKViZ51/Pc5lFLPG2r1qePx4e9cfB2yJBaXOP854pjY1KEB9vHJ3BIUbL1yT3ZoHizAMZ5jcMA4DpHOeLCRIevYMImPtemThCXG8QghYFcf93bp/SIpqyso2H18s5+/Zyba7wYErzW33Qpq3InVv2l/3yL6Op9dEfnMPJqIul4SmweGFyQ+cxsHB+cnTrPPZfpwqwsDn5G/Ddw0DA4nnxPf0jK+cn/4wKuvvAwjkQm3ipxkEfDM0FSyGqAWBmBdA43cBgAWh0AUaTld5dBVqngFEBqjYXPmoApgY0bDHDqAfUWta9jS5Aprve7irQbFv1VnR3uaduMmpWP5q2vFUzrKTR1cakHs2jz+iuAo0OGDGQt91GF6iM/PfXu0DFsPrqKu8B1aH//lofqA385dUBUO8nL68MJ+PJlwPAgtvngAYtwKgDjY57PMd1jaz3FtVV71nXeUruXGudusajMmxC1xvBfat3J/PCe032Tci0v3t10v/uCgAK6B8DFtzKfOotgzEt+H1l8yXsnURzw+k3v25k7RhVYLAoHx9uPpl8uaYDzY7//v/tZWSFueb0a06MdYHHXxYwUjpZmDahT0LwLRD63RWYhrO3e5RPUhHJ5P8AUPUQv+R8yqTuipdAZwXLATl3FPms6e0bVdzEUh+pW4DlAaJ8TzxjKDkfzg7mltPHr/4aSOuhvvwsUK8AGwvANsdBXV0F7h9PCMHVVf/1kQE4DhxX7R3ee49hArZXHbu0DFS08T2HzSpWNxfdz6jg/CKwEz2BhAxDxlDp74NV6qif3JlcWHoCaN8FqzQn7vSbN4H9t4D6MsCbQdr3AgAW6kB3aN17csfKB8u742/eBA7fA7QaMDpxlZOobud+ADB6QLMK9HWr/Oh9QPeMxeZN4Pi2JYLhrzn9GbbhUlhu3gROPrT6yFsFbd4ETrfdfXHKO4+A3o6/vLsDdB9FKN8Fug/95f1DoH3PXw4Aez9yFbPlJ4DGGnB8BxidwrV7bd4EBkfA6V1xXe37wODYrWB3ru2/5da11FeBlSetE2r7LtA6PxkD55nTu1Z7orYk34SMIXD4NlBbmcyHjRuAVpX3zzSAA4/9x8pTQH0Z5PRP1IfDd9zf2LkWNN6nd63Tq+mZG8bAqq+2BDie3eN23nXP+doysPoUcHoPGBz625HNp1EXOL4l7ldGmGtOH0SoaBQ7T2tUeFuJa7OdusiPCEZrE2Z9BaOlq1xDFmfnMqUbG4XLuZts0k1O6vSGIjari96bnRuF1UTxUYiOfBS5RnOiiVXKjpYCyBsCIWXrKLMqS4dWdESR6RdfITKfRH/YAQ2PQQBqFcrES9IB/4m9RF4r4KnQsJOqA0CtY/nM6XwaQCay2YSkzCmO+KIyYiIoJzBhQK70dsZwxfL0wc+riH2KkHJvEiY5G8LFKhOvaU1qnSWHahhnkMj5L+CdlKyCgqHcNw+ytzbzYz6Jfmdv7FIexOnXZWkGBVCJm294YmvEPWHITgiNWjqfq795E/rCJfsvAtMmi0SThIh1Qem1ZAup+JxQESDW5+Q9djE2GNXTgbOeImxKSmEiABi1Bfm1Zhxb1qhMTTb29WlhbmX6zvyrapo00FElhBVPykzGjeIpyxrUrFZSCSEBAPrCFvTWOWshcZx0/fh9e9KKx8asrQCwrSBSlkMxaaR82XdSX4zRTUvTQkRxkGBMXVwkke/R+CaWMR5M9M2nf3Iyqwu+JOhSiDzW9T5YtYm0N2CmNcaB3rLGfHL6k1gCY8KukkatVo0+Kfl5kRYnLm0rbU6PNNeRe9KQ4ec0bUuo0cJ5Lipn3P5IuFihjD5rJyymFEr3bCItIi0gnqGB0IK+b5CoBtl6aiv5TERnUozGaqzuxMF8En0bBBp74opEJgXyv5k6+ps3MVq8MpaHsopX6eYoeDXodrQsK0lGhCM9X4/oin0sVs0sJq0/qD+U8HArzVJulaenx0mP4CY3/U0ffHIb9f4lWLABiz3sFMi0COKaNIhKxoRpTon+xM7dkdurEJKgsc4qoYEI0/LYM1rn0D/3cfQ3PgGjue6STToxcphWGS+9NLNP8XC/fpCSLToxGy5dVu1FrPpzMhTLGNZLRN0soojQ4obCngbMGhfBNoYCOhBTWOtzSvQnp7yKRiDCmDuNVsfkg8jk7K77U/p+UTJqZQKbs+mf+8SkyMlkxXHK2jDI/4Apuq+Pbx9z6CJdi89sM+pg21ZGPmITy3JjWpxz2noJNo7KmjdkUVQDn5FkwcoavG/WOCdwkIQpoTNonLGJgrlV5FqwvlatIlfmJgW/MRz3RqlExiuS2Km/eRNgJiqDI0ueyXF/leExzOoiGGl+UiSTrSqElmCev+QK3uwwjkMf6VuI3H7VFbaB9ztXtRrI8FqxCJ4xjfHmHYiIlibJdCnuuTMSlAdBk6TQdCNEVJPHRJJaXlmWcsKw5/zj5hAM2fk0zCWn70xLZ+irFVLi1JNClF4vTUxtMyANRnMDhh2ueWByJp56x7Iu8cjKrSBh8ullWQFJroleVLSQ4izgiCILqZI3qW5ABAVOVp2L5MMQB+k4AgikT69jl0cVzQRO3DibX8oLIYKvSZgCOkzEZTTWFPpTyvTjw57QNU1LnAeWT2coQ4EY9Ewx2HgOg/WPQW9Zgeb01haGyxNP3+bea4DMqcqGOMGKwNJKao0hJxZRxEphi3Sw9hFBu4TB6tPC2pIhZAZlnFrQ1Q4Asyrb8ApINkI3zHROKLwDWyzRmO/75U815lK8490oqxXCqJ9MJq+SiCVJCkXVfhQFrFKHvngR+uJFu8CTnIPp/ty6Y1AhJn8w7M1Alv0rSf+L8oHT7keCTUlfuJAo/wOZo2hWNgBgjqz4UFKEjE+Ezc+ftlKOyuAIZoYmnAXcslOAM7j2v9UKxZrfSpmzMl6/2qx8IdKsdI3Nc+NYPkZz3dIJ2Gju/QiArR9Q5optWagvno6Eox8vLu/A5cQli/qjiPiWIWHEKWa1HLRRh2sqXoVCeXoKIRB4UBQP35zzFKhC1fM4LuaT07f/ddR/tZmhnH7wNKMg/GEgRktXQu8hcwiz2nR5RpI5hGbIp6Nqvl0CA4OG0eKlcWwh7x3uigVl0nsZx7GJdAwexW9kCyM2eSatDSotrmRKpxO9dR5Vb3RKG6y6EJBhLQqsb+qKNquENHerSV2OsjwrFmV2qaEKOE5/VlEUSUAS9Ddvor9xA8PlJwAAw9WPuJSBTTuUsG8Bs8n2nTuEAy+3xBEm0ZbkC7YIS/QPG8UOvjJUDQAWVN/01o0znmKPcf/YhZ2SZEnOeUgZFlEgQKmTHiD1OA90AJucZNPZyOSYS6Lv/R61itpr8uaXXhn+PBDfqUKrwmysob95E2ZtEWZ9SXgbH/q5cXTLMm9jBoQcujRsQ9bEKqo8331/HmKFMLNAB6qB3TSpfsbzXNTYZGNiGCMEikKu5eAKPGvca76aYNHrCwrZ8kwDzlibos0tI8yleMeBE5e9XtGsk7fvG8o/6sgwXVdP+9nK2VSgEulzljCW9zMGbXhseT6Shlp7GwBAhkVoKsNj6NXW5EFTBzHdTzwZg5zwyKx6ZLJ+gT+BU78Uit6onN+B0TyHSn9f6Tl5hfyJiHnKgvqlvkFqegfAFqamGxGBV9xOY23EGgf+mems59w4fSJaJKLfIqLfJKK/m2Vb3qHUNEJDEEZZ9QRQImMQwWysWYuYCP3Nm2CV1thO3GisQefyADQP3gSYAW3ozqHqWBDJFbyK31uBa3SFox56slmFQDOGE/GF1OY92Smm2j+waykSo5DuZqG3tsa/q709+1fwpsQnDJqcDsLMfKOYDEuQywlUDYmoHhF9lYh2iOh1T/mLRPQOEd0ioi/bxX8TwO8yxn4JwN9I0q5Cx+x/Jq8n8pTly8z5CJoyNxisP4PB+rNWILjlJ6Smk0ZzkgCmefAGyOj7rB+0Yds+FYg4dz9k8tqJe7zb5FRmbUHjk4FHvGOIRSUWoU5L8RpwMlXmirl+J+KkVZ5VVKa7auWvKYrMOJpQ8TINvpbjvHPMk1COp6ekrO7XALzIFxBRBcBXAHwGwA0ALxHRDQBXAThJKjMNoiGan8uNIHtcoDsysNqa3KMSilnWVlY469tSf/Pm+L/B2jPon/ukUPmm6R3XIqqf3IYj9uFtuSuDYwSKfTzlleGJ8INXe7v2L8/9/QPfvYAj07c5fa6f7gB201ei8vCGQDAa6+k0Fqj4DHwoTkOeIo9MP9M0mxLMWsA1xth3AHhn9qcB3GKM3WaMDQF8A8BnAWzDIvyB7RLRF4joVSJ6dXd3V3abEvhcq4uNcMcNPpPWMIewDSXig1VbY67Nsg76BIYrT4NpDYyWrvmIEhl9SxzEiWVqXHLtwepHx78r3sTWk1aFpcQMiUXJhIgPl5/kqpHI29lI2kZsCIkq8/wbrUIGgs6dsLzXckNi7tizxqX1RbfGEV0rSojrLITaVzDh6AGL2F8B8HsA/hYR/UsA/6fsYcbYy4yxFxhjL2xtbcluC0b7kfUvd+yueuT33ry53s9BoLlTnM41tArM+hIGG8/ZYaHJthRaGofGdcJGuB4btQGtAsal2Kt2HwudhiyOXqTgtU6F3nJtdDK26/d7WAbPLalVjURJPVWoRr8MJHqTd7DyNPCIQqYimLUOjuxHJCKZBOuf9ysZ550oCHKz3mGMdQD8Ql7tAZwXoQBeEX5PUZzjRerZrILaKs68mRkMBTFy+ps3x74BZm0FuiAIVmV4DFb151qtDA4lYgDr44yWrqHWvucpF20S9onDQ2zIZYmjAI6xMWuLagm6mamu1E4DYQHXRFnauDKz2rKthzyI9Q5Zc9vWaccUzB1fT1y6oNmR6YtwH8A17u+rdlluUCaO3H29oVGkzbhExnB0A8PV65blEID+uU+NrxvNcxi1xCdNkzcftVEZHAJUsU4ZPJjpi4xZGR6ntoPzhCPIa5l/NychTnjlquEMvErjNAmYKE6TBRYp0mm08XadNmKGCB/fyqVv5WGdNPInOllw+t8D8AwRPQWL2H8OwN/JoJ3IqFe1hKn4wpGtSKjcleKgViWMdIWxs0VCXvD+BOCih/Y3b6Jx+B7I6MFobLjMAR3I9AMOsTalaSklvgOBYXvXJyILvj5yE87K8Ni/OYn6ILzk58xJ71unH9thrr9xA82DNwLqV8C0j7WmjljrTTQ+nOK+CGabSU02vw7guwCeI6JtIvo8Y0wH8CUA3wbwFoBvMsYSzoB08PTWYvhNNo57w5LEzglSS6BD5BMpDNafsc1Kr8JorgGArUuw8gWYtWWMlq55a0JlcACQW5cAAJW+s0l4liZjY6LKZxHT9N5YJGTW1Oa3NupNiGoqykXP+LqcphSi23oIod/7VyYCk9QndKqTQxTeIk3/BpmYWRZTKGsk4vQZYy9Jyl8B8EqSupNg5KRG9ARa452xwqxzekPxdZ5rZIxNnSEpUUwMV6/7yvqbN6ENT6CNOmBabRxXZrR4ZRwcziEErOI2MW7uv2b/WsBw5fqYk66173Eu/6oE3D23XXoI1RARvKmpMZJayWujDhfeQM0008rIxovQImxMYTGTlMMbH8OoL9tPq1vvOAH/XGVTSvMow1yGYdjvWh+ZAuJZdAfJwyroOTt0NaoV6WZUYjZg1ldg1t1Zw4zWORitcwAzx7JkVm1iuPIU6icfWDdRBWCGRQw9cePHpwMid+RJZgp1CmAGR9w1GM31MdG3vFSjiSKq3YfQF+SWdpG5Zg/hHi5f84iLnFOKXFBhNDdR6e9N6osc8XQSF0f9GVMo3qn0D2BWm0WQ7ACY04BrIyeLU0BSBW8Ulf7IUPrEYXMnS84/702mRM4gDazaBKtacn6zvjxxSDv3CfTPfWqcsrK/cWP8GKvUYTSsMBWDlevj8ub+67ajmn+ZV+14P94Noda+j3HIAEFegnF6TCEkDm0RwTuCMSLBOg4XTRncxqoWLM4NVWW3EzkWsHQbQgT0k6awpOeS6MfJRt8fReeg8xbtqHoJl5hT8MRDq3IWSE9zjE51TIj0hQswGhsTRyqaiIyqXcuXxZttShu1bRGT1ZYj4gDsDUEGxiyFpddSqb8/OVX4ThwM2qgLzej7nb3Gz4THIOI96ccbD6+49gS08yncRXFxeIukoEikXGROclhJz+3RYvRPsFDPRhAzl+IdZ3JRzW9a58AwGSpa8HkrjKh3h7rP6atE8SCOsDq/MBtr6At8D/rnPm79YKYlxjENy7MZwGDtWTSO3rWer61Y4ggA+uIly8QUE+JpNDd95qGOzkEb+WPBOyIqr2x8rKdg8HHzk/j4npPGsA33GZ2w2qpZOazH1lUexbDRh1FpwKHGZnVB0Z8hYswg20qHgXx5AKLkQgCAK+utUPoUF3NJ9M1KC/rKdWBFHtN6vz3E+ZVksTbyJiRniXCVyBCkWSaWvJFNtSk0V2WV+qTc1EGmbhE00jBc/Qhq7fsgU4feOgdt2B7HQhotPeEKc+HUNVq8jFrnQWgXnXu8p3YrjtIEVntWJJjmwdt2qeaKjU/GwDq92GW8c502OAIZPWESeEc05YvVb59CNL07yRENoNZ+MKHrglONUL8y7mSczF3xMJdEH7AUZpRBmsRQmX7qLfJ1l1Q/DsxS950OtOqYw19frOGws4jB+rOT61ziEKO5NjZh5WG0Nsd6CZgGyByOw0v3N2+i0j8EGQOY1SaIYWzS6pxEzNqKFe2UnHY2ANMi+malDkINRm1xrBfhYY5Dda+OFdf1U2tjYoIE6WMfC8ZceQ+a+05QYY+Zp2lFYa0YA3iDc2ujU5DTNp/GcxzCQwPT6iXRTwoKOU55CXhUTrqvG7mmYiyJVwkAaNW1qVtx1dMQa2oVMM0tgpU5jclOIgDALv0Y+ppfXOO6n1/cpKF/7lOW/TwRYOpjHwc+RIdRX0VleAqzsQqzserSDTCtgdHiJQAWg+nL7SDAxlId97Uq9IULVnwnYLIBUMU6Qdm7Re3xD4EnbgKCE0hSzDXRT4wULHVMk6FR0zCIoSieFpaaVbQLkCmshBibSw3cO4hukXKm4VXGEklTdvKbxUhSzmPoWEyZOsgYjk8Z/c2bqHZ3LJPN2gLY6hYwgIvom7Ul66RDGkaLl1Dp72FjsY6lRlXdZyIi5pLoO2KQLCKZ8iKWUthSosT84dxSHftttRzDLnDiLwf6wnnoC1Z0V9ZqAgNLfCPcQOwwIOeueiOypou5JPp5QRRnR1yWR2/SQ0F8SEqUCMSFlQZGJkOzJvbHadUrViDFgqAo62qu7Q2zzllw2BmFcvvGrFH8HLG2EJzNrERx4F1Lawt1VDTChZXG1L5jo1rBlTW5WXYqeoc5xFyOyjiWlGBv5RWvXmuYLKxjBhlH9ZxVPLGxgGsb4THHS0wfT24u+Myb61UNNy6v4PxKE5tLU0gzCEAvLRtiYS6JfhAurvhNuWQQbQLVDMxAi4Zcsrql0MbFVfVvOU8Is0pLG8uNYkqBZy0qCeWysMIx1xRMNMZul+3k8Jl9plr72UJUMUFB1lDuWGnlT4S1gg/2UrPqDaobG6qvmpbH7IWETqJRMddEXwQtI9dmGbpDfeacqvLgJEULq4gy/iIe7KbBMZ5brKfOMKWJpzYXcxcXpvYZcv6cBZzSyZGl7jSMgHvbniX7/GmjiNzkUoqijYWG38pE9spJnf4WBW3FBRGBiBKHLckbUafTlfUWnrkgtt1Po/6izO65JPrLTWuhhhGRMI/c0vAmP3zq6ioW6tEIVa2IbHhBsLYgzyWRNmope6VrGiIR37SwUK9IzT/j4imFbH1562jmctVcXW/huYvL4VE002jMJ9Of/Z1iWgx3VLGFs7kXAUXqS95IO9JsRSMl4uudLs2q/JmNpfBNMIt5n+ZJMS3MJdEnItSr4a922veGRUoOQ2BSUJ4Y5gsiMcf1TfX8y0VG0TevTwV4qwateRWbfUcyUCmgmDFNzCXRD8OlNcvUz8l1myaGpV3+XEG0YVcVjQHWFmqu02YRuT4vgnQJaYs+wpCU9saVufPGHionhKz6kxXOJNE/txjvQ85KIvSkJmCi08q8IYmidF1RXr613ECL01OcX54tRei0EUYkw75gVJ8aR7zIb+pBHr8ixFEz5S0SzpXoE9HPEtFvEtHvENFP5tm2px9C8zPv0CvlzPV69RZAGbya0PQxqafjNDka1Q0vSR+jmP3ysZiK4pwThCL3sVHLh1wlUYJvBWzseStsZVAeRSL6KhHtENHrnvIXiegdIrpFRF8OqoMx9vuMsV8C8EUAfztel9NBHJvjOAT8uDeaidPBLEI4rARcWJUvvOcuLkuvFQlX1qNxmFmgUdNweS261/OzF9O3vDm/0sCzF7L9dlFJ8qyu6yhb59cAvMgXEFEFwFcAfAbADQAvEdENIvoUEf1fnv/Oc4/+E/u5qUHk1SiKkOm/J/h6e1DGoQ8DEbCQoXx4HkJltHKWn4vQqGo4FyOuTiPAisbBvIbQWAxKZl4MRl+d6DPGvgMnGeUEnwZwizF2mzE2BPANAJ9ljL3GGPsZz387ZOFXAfwBY+z7onaI6AtE9CoRvbq7uxv3vUIhcjaJI8v2OWMVWJF7bWP63CMAfPLKaqiZn6ZZzkxes9sweh7nCM2PS5pewQWWlEhhmizRPFE5GTy9tRgoBikaopx2ZuGbJ2WJrgC4x/29bZfJ8A8A/FUAP0dEXxTdwBh7mTH2AmPsha2trYTdC8aFlabLc7EzcMfe9nL+nRkMqcBDxYy1KGhUNXxkawnNiHLcpIvu6noLNy6vJKrj2kYLV9dbShzvNBC0sR11R4lkz3FOBl447W+lUFcacz4vq6W1hVou4sdcbcgYY78O4NfzbDMMV9ZbePeRlV8zjNNXMccc6P6kDUXdKLaWG1hpVfH+TsdzJT92hUi2MN19eGJjAdUK4c5+R3KHvCxafwhJHUw1jbDSLG6cmlngRgHg/EoTj08Giep49sISGAP2OzEyYaWMsHGvVtT8i5IiKdG/D+Aa9/dVu2xmEMSNiUh1mExfZvhSq1ImfgFJcGGlIbTWUNFtpIVPXlFLDVepEBY9du5xguc55nibSw00axVUNUJ/VJzsSg6ypMtFttDxYm2hhqNufCdKK25Qsj44q6FR07C+UMej436yCqeMpNvK9wA8Q0RPEVEdwOcAfCt5t/LFJyTHeZWwycUi4xaSmob1EwSJO7/SwFMpeqdGHd8wgqZphE9dXcXmUgNLjWrmR/e4kSmLOK+iIC3diMP5Jt2m0jhtVzRKNZDdtBDFZPPrAL4L4Dki2iaizzPGdABfAvBtAG8B+CZj7I1supodNI3w5KY/LKvpofrev0VoRQwalgTnVxq4Luj3NHFhpTkmpCKlq6qttUO7Hfd5xzXeu6HVqtlwrdc3F3B9cyEWqeAtb7ynkyIgDz6/NoOpCqOeCDZiOnlOG8ozkjH2kqT8FQCvpNajKWGlWcNio4JaRRsfJw88csDHx4PQiSFLxFzVCOEZdaNhuVnFQpCJ2JRQ0QhPnFtAq1bBO49OAQDnlurYb6vLVZ2Q1FfWWlhp1nLdTAFg2ZbJ73F93lpuKOU8nhWCt7FUx0GEbxKEC6sNHCcQw0RBRSMYJgs1zeXXahYRWS+sNLEToHPwnnaKIlQrHsWYIp7espxKjrrHAADdCA+epmrbv9ioojecvjIpL6y2ajA5xfjGYjSi7yjVNY0SexinhTDb8o+eX8LpIBrhu7a+gJ3Tvk+klq1MP/6zV9ZbGBl+8d/55SbOL/vHJwtd1scvLeOkpxdmXojQqlcKmwN6NliSOcEM6c+EuCSxVw6KfOhF0iHIwrrBsbSJ6wXriJha9YqQ8AVhdaGGZzL0NA2ac3G+xcZiHRci5JnmRSDe05pKUhYnx8ICJyYjyp4ROL/SwBMxiLZz+ihytNKS6M8w8jCy4RVXmxHtprPY5NKw3fbCUe7Oqow2KoIU/Us5EiuVaKXLzRpuXF7JPULphZXmeGNx9EoqSvl6VcPHLi1H2hjzRkn0BcgidkgWUJEvJ0XehNBxxpIusBk/LRUBzhiLQj3wnrJJrMDS9E2Jm4A8KKlKFNSrGj5xeUWZ6QnS6RThtF8SfQEa1Uqg4pAXMahO7Swi7BVg/qQOh9OUjX9jhryKg/CJyyv4+CWxWCdrs8ClRhXPXFjCesiGnibhFsW6kiFpfmDASreYpgjI8QmJo6R3LNaK4h8xHysoA3z0/JJ08sVkPM4s0pzrRST6cQ5cmkaJ0gwmkRkTqaUjTAufvLKCJ8+p+26oJqkJgvf9REnpvVAhynGIftyTSlYo3goqED5+aSVUSaliu58mkiz2aYaCbVQ1rLZquLpeTIuGeULS5COT+9IhVkXgcIsaB2kaKIm+Ap654Jbxu0L0RCCkT55LTvB4rkG2mAqwxnwgsm33Q+ztnY0pjXdw2koiLkgrxHHUwHFx8bFLy/gYF7QrrkfwNFGUZCNpI+yt8nrvkugrwHtU5AOvnfTD4+czZhGx5YhBuDaX67Hku2FH0HpVw/XNhfGmoXL0nUU4GySR5RwWB0mJppNsOy9xSq2iucRGRWQASkzA+4zlFZixJPqK+NTV1XHI3Y+en3D+jgduFvHBL622xg5jaaKikWsDenpz0Rd/6OmtxcS5duPAIZJaRGoVJje9vNaK5E8QBUFxZupVDR85vxgp16rMYmptoab8TYKcBud5I3j24hI+JlGQO3BO7pfWmrmkYHTmsnfcn7u4jI9dTBbGOw6K60FQQFRse24RVBfS5nIde6fRPHOJJmKPLOTyokiEi40qFhtVHHXV+vrcxWVpbt0oBPz8cgMaAesRLS9a9QraCqeutHHj8krohhMlVAY/v9YXa1isV7F92AOA2B6e51cageEC5gkqsvtaRRuP8+ZSA69tH2fap2vrLRx2R7554HU0LMU7BYdXzp8l9ySTSQc12arnaytcr2q+SU1EuLTWxNNb6pYbmkY4v9JUVv6JNkGHBldyYGmztMy4uh6uA5GBH78LK03X6ZTH5bVmphnVRITM+1nWF6erd7i81sz0O1YrWqAkoFGxvnFeSY5Koh8TzVrFCtG7bB3Fw5JmBE0pWXgDEU76arFdLq5aC3na2bKcuPVpQWUzWFuo49JaU8nNf17hFe/INo9zSw2sLUzHE5liivLiQhZz7dxSI3G2tCRYXajh6a3F3BwhS/FOQlxabeHSais08NpQEKTKweZSAw+P1BIz5GF2OQd5xbG51FDKdFZ0qG6Yz1xYmjlZ/fpCDQPdwPnlZqRgfHExTZPlMOQZgnsOlncxQESBIpXDjsWhR41fA6g7qzy9tRgpibMMs2gyt2Zb2SykHII5bVPLJ84tpCrOqGiWPqZZqyjboudhN++IS4KmLhHh0mqrcM5LWaIIgdim34M5wkfPL4MxhpOejs5QF3IvcTz6njy3iLcfnvrKZcrXs4jFRjUT65y0CeRqqxbZDFS0SThM69WNllC06MwD1bYWGxV0BumljdxaakAjcokszq808OFe90w7Sq216jjp5W9wwONsUogM4YR9XWlVMdDNsUVJElt42UaRVQKVWRMTALPZZy9kYrU4XsyOzmlaICLfqXalWZtqn7KETFFeRJREPyMQ0ThXbH9kCGWzW8sN7J4Wx5RuuVnFaV8XhoUtwrFUhFpFg26IOdS8nF3SgCxPswwXV5q4Z3SxWMDMaWcRylZWBWBOyhmTA7wE/8JKA8t2CkBvNqY4dvxPbS2ilkJkQodbFslYiyp3vX5uAe2BPjMpCmXQIo5vq17BsykmX5md7TE+ahUNQ92cLt0twECXRH8KOB+QYOHSaisy0RclmHCUvwv1ytiKxbH3j8IdFtniAbBsoKdhcnhhtSHNh1yimHhqcxHtgR55g503lES/IFhbqOG4Z1n4XNtoSU31VBOMN2sVPLm5gKV6dZzovVGt4JkLS4UMTzxriJoWscT0Ua9q2KiejexoQch19RPRIhG9SkQ/k2e7s4BrGwv45BVLybW2UJcS/fUIXO1Ks+bjapq1itQixUmqwUeWLKqCNK+olWcF4/gwU+7HWcE0gxwqrRwi+ioR7RDR657yF4noHSK6RURfVqjqHwP4ZpyOlrBQrVBm4XId64ppe/GG4emtxbGSXAWz6HeQN66ut3B+pWGZetrDVdQNfx6QRqKY2G0r3vc1AP8CwG87BURUAfAVAH8NwDaA7xHRtwBUAPwzz/O/CODPAXgTQHkujoGt5QYOu0PUKhqubbRw2cxnGMNk+lmn9hO3GU0qOUtWPFniynoL+22xtVitoo2ttlaaVWwtN7AZMxx1CTkcK59phb4AFIk+Y+w7RHTdU/xpALcYY7cBgIi+AeCzjLF/BsAnviGinwCwCOAGgB4RvcIY8/nJE9EXAHwBAJ544gnlF5l3XFxtji19iChSYhArln866pulZnXse/CJyyslNzhD2FisK8V3ISKfVVmJdFCvalP3VUhCCa4AuMf9vQ3gL8huZoz9CgAQ0c8D2BMRfPu+lwG8DAAvvPBCyaKlAEdXkAZ409BZsYIoxTslSkyQu/UOY+xrebdZokSJ2cFSQR0B5wVJRvc+gGvc31ftshJzBEcGGRY6ukSJNPDxS8uudI+zgFnbpJL09nsAniGip2AR+88B+Dup9KpEYdCsVfDJKyu5RGZMG45YZ8ZoyJmEk8N41gj+xy8t55YPIC2ommx+HcB3ATxHRNtE9HnGmA7gSwC+DeAtAN9kjL2RXVdLTAs8wc8rwXcaqFc1XFpr4slz6uadJaaDy2stXI6QR7goqFa0mdFtOVC13nlJUv4KgFdS7VGJQiNKwpciIE7+ghIl5hmzdZYqURhk5SBWosQsYNp5fZNgtjQQJQqB0j6/xFnGtO3sk6Ik+iUiY9ZkmEnx9NbiXOTbLVECKIl+iRKhsNJQTrsXJUqkg1KmX6JEiRJnCCXRL1GiRIkzhJLolyhRosQZQkn0S5QoUeIMoST6JUqUKHGGUBL9EiVKlDhDKIl+iRIlSpwhlES/RIkSJc4QiIUlQZ0iiGgXwIcxH98EsJdid9JC2a9oKPsVDWW/omFe+/UkY2xLdKHQRD8JiOhVxtgL0+6HF2W/oqHsVzSU/YqGs9ivUrxTokSJEmcIJdEvUaJEiTOEeSb6L0+7AxKU/YqGsl/RUPYrGs5cv+ZWpl+iRIkSJfyYZ06/RIkSJUp4UBL9EiVKlDhDmEuiT0QvEtE7RHSLiL6ccVvXiOj/I6I3iegNIvpv7fJ/SkT3iegH9n8/zT3z39t9e4eIfiqrfhPRHSJ6zW7/Vbtsg4j+kIjes/9dt8uJiH7dbvtHRPTjXD1/z77/PSL6ewn79Bw3Jj8gohMi+ofTGi8i+ioR7RDR61xZamNERH/e/ga37GeV0o5J+vXPiehtu+1/S0Rrdvl1IupxY/cbYe3L3jFmv1L7dkT0FBH9sV3+O0RUT9Cv3+H6dIeIfjCF8ZLRh+nNMcbYXP0HoALgfQBPA6gD+CGAGxm2dwnAj9u/lwG8C+AGgH8K4JcF99+w+9QA8JTd10oW/QZwB8Cmp+zXAHzZ/v1lAL9q//5pAH8AgAD8ZwD+2C7fAHDb/nfd/r2e4rd6BODJaY0XgL8M4McBvJ7FGAH4E/tesp/9TIJ+/SSAqv37V7l+Xefv89QjbF/2jjH7ldq3A/BNAJ+zf/8GgL8ft1+e6/8zgP9hCuMlow9Tm2PzyOl/GsAtxthtxtgQwDcAfDarxhhjDxlj37d/nwJ4C8CVgEc+C+AbjLEBY+wDALfsPufV788C+C37928B+Fmu/LeZhT8CsEZElwD8FIA/ZIwdMMYOAfwhgBdT6st/DuB9xliQ13Wm48UY+w6AA0GbicfIvrbCGPsjZq3O3+bqitwvxti/Z4zp9p9/BOBqUB0h7cveMXK/AhDp29kc6l8B8Ltp9suu978E8PWgOjIaLxl9mNocm0eifwXAPe7vbQQT4dRARNcBPA/gj+2iL9lHtK9yx0FZ/7LoNwPw74noT4noC3bZBcbYQ/v3IwAXptAvB5+DeyFOe7wcpDVGV+zfWfTxF2FxdQ6eIqI/I6L/SER/ieuvrH3ZO8ZFGt/uHIAjbmNLa7z+EoDHjLH3uLLcx8tDH6Y2x+aR6E8FRLQE4P8A8A8ZYycA/iWAjwD4MQAPYR0v88ZfZIz9OIDPAPiviegv8xdtzmAqNru2rPZvAPg3dlERxsuHaY6RDET0KwB0AP/aLnoI4AnG2PMA/hGA/52IVlTrS+EdC/ntOLwEN3OR+3gJ6EOi+pJgHon+fQDXuL+v2mWZgYhqsD7ov2aM/R4AMMYeM8YMxpgJ4DdhHWmD+pd6vxlj9+1/dwD8W7sPj+0joXOc3cm7XzY+A+D7jLHHdh+nPl4c0hqj+3CLYBL3kYh+HsDPAPi7NrGALT7Zt3//KSx5+bMh7cveMTJS/Hb7sMQZVUF/Y8Gu628C+B2uv7mOl4g+BNSX/RxTUUbM0n8AqrCUHE9hoiT6RIbtESw52v/iKb/E/f7vYMk2AeATcCu3bsNSbKXabwCLAJa53/8Jliz+n8OtQPo1+/dfh1uB9CdsokD6AJbyaN3+vZHCuH0DwC8UYbzgUeylOUbwK9l+OkG/XgTwJoAtz31bACr276dhLfrA9mXvGLNfqX07WCc/XpH7X8XtFzdm/3Fa4wU5fZjaHMuEEE77P1ga8Hdh7eC/knFbfxHW0exHAH5g//fTAP4VgNfs8m95Fsav2H17B5ymPc1+25P5h/Z/bzj1wZKb/j8A3gPwf3MThwB8xW77NQAvcHX9Iiwl3C1whDpB3xZhcXWrXNlUxgvWsf8hgBEseejn0xwjAC8AeN1+5l/A9oKP2a9bsOS6zjz7Dfvev2V/4x8A+D6A/yKsfdk7xuxXat/Onrd/Yr/rvwHQiNsvu/xrAL7ouTfP8ZLRh6nNsTIMQ4kSJUqcIcyjTL9EiRIlSkhQEv0SJUqUOEMoiX6JEiVKnCGURL9EiRIlzhBKol+iRIkSZwgl0S9RokSJM4SS6JcoUaLEGcL/D3dQxci3eknAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(s1, alpha = 0.2)\n",
    "plt.semilogy(s2, alpha = 0.2)\n",
    "plt.hlines(np.mean(s1), 0, 20000, colors='r')\n",
    "plt.hlines(np.mean(s2), 0, 20000, colors='r')\n",
    "# plt.yscale('symlog')\n",
    "\n",
    "print(np.mean(s1))\n",
    "print(np.mean(s2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
